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1 .  Introduction 

Despite  significant  advances  in  primary  and  metastatic  lesion  detection  with  static  PET,  the  ability  to 
accurately  detect  axillary  metastases  at  an  early  stage  remains  the  greatest  challenge  in  breast  cancer 
imaging.  At  an  early  stage,  radiotracer  uptake  of  metastatic  tissue  is  often  weak,  due  to  the  relatively 
fewer  cells  involved  with  disease.  Moreover,  the  diseased  tissue  and  abnormal  uptake  is,  in  most 
circumstances,  embedded  in  severe  background  interference  and  count  noise,  thus  hardly 
differentiated  from  surrounding  normal  tissues  with  visual  inspection  in  either  2-D  or  3-D  static  PET 
images.  Once  metastases  develop  into  fully  visible  lesions,  detection  may  already  been  delayed  for  a 
certain  amount  of  time.  To  avoid  or  reduce  such  delay  it  is  desired  that  the  molecular  changes  of 
metastases  in  PET  images  can  be  detected  from  interference  at  an  earlier  stage  before  they  become 
completely  visible. 

This  proposed  project  uses  a  different  approach  to  enhance  the  signal-to-noise  ratio  (SNR)  in  PET 
images  for  early  metastasis  detection.  We  have  been  developing  an  intelligent  detection  of  early 
metastasized  molecular  feature  (IDEMMF)  system  for  evaluating  axillary  metastatic  lesions  at  an 
early  stage  with  dynamic  positron  emission  tomography  (PET).  The  objective  is  to  extend  lesion 
visibility  from  limited  2-D  or  at  most  3-D  to  a  higher  dimensional  space,  where  molecular  features  in 
the  malignancy  and  normal  tissue  can  be  better  distinguished.  Thus,  increasing  the  SNR  of  metastatic 
lesions  not  routinely  visualized  becomes  approachable.  The  premise  of  this  approach  is  that  molecular 
features  available  in  dynamic  PET  data  in  the  primary  tumor  of  an  individual  patient  can  be  used  to 
identify  early  metastases  that  are  obscured  by  overlying  background  tissue  activities.  IDEMMF 
“learns”  the  molecular  features  from  the  primary  tumor  and  normal  tissues  of  an  individual  patient  in 
order  to  identify  the  tumor-like  multidimensional  (temporal  +  3D  spatial)  molecular  features  of  tissues 
that  are  undergoing  metastases  but  too  weak  to  be  visualized  due  to  server  background  interference 
and  count  noise.  Using  IDEMMF,  these  features  can  be  enhanced  via  a  combination  of  sophisticated 
PET  molecular  imaging  models  and  advancing  image  processing  techniques. 


2.  Body 

The  whole  proposed  study  consists  of  four  tasks:  Task  T.  Developing  the  mathematical  formula  to 
linearly  map  and  identify  the  physiological  features  contained  in  PET  dynamic  sinogram  sequence 
(Month  1-8),  Task  2:  Developing  the  schemes  for  objective  reduction  of  dynamic  sinogram  data 
guided  by  the  identified  TAC  subspaces  of  the  desired  signal  (tumor)  and  the  interference  (normal 
tissue  background  plus  noise)  (Month  4  -  12),  Task  3:  Deriving  and  analyzing  statistical  hypothesis 
test  criteria  to  test  the  presence  of  an  axillary  metastasis  in  the  dynamic  images  reconstructed  from  the 
compressed  sinogram  data  (Month  13  -  24),  and  Task  4:  Clinical  Evaluation  (Month  13  -  36).  For 
each  task,  several  subtasks  were  defined  (see  the  SOW  in  the  grant  application  for  details). 

The  hypothesis  was  tested  in  the  project  is 

Hypothesis:  The  molecular  mechanisms  and  kinetics  of  [18F]  FDG  in  the  metastases  of  the 
individual  patient  are  similar  to  those  in  primary  tumors,  but  differ  from  those  in  normal  tissue.  The 
multidimensional  molecular  features  extracted  in  visible  primary  lesion  and  normal  tissues  of  a 
patient  can  be  utilized  to  enhance  the  tumor-like  features  in  faint  or  invisible  metastases  by  cleaning 
up  the  normal  background  and  statistical  noise  superimposed  on  them. 
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The  Intelligent  Detection  of  Early  Metastasized  Molecular  Feature  (IDEMMF)  system  under  the 
development  consists  of  four  major  modules  as  depicted  in  Fig  1  (1)  tracer  molecular  modeling;  (2) 
primary  tumor  feature  extraction;  (3)  feature-guided  filter  design;  and  (4)  detection  criteria  design. 
The  functions  and  interactive  relationship  of  these  modules  are  described  as  follows.  In  general,  the 
tracer  model  is  the  vehicle  to  define  the  relationship  between  the  noninvasive  in  vivo  measured 
dynamic  molecular  PET  data  and  a  set  of  parameters  reflecting  local  rates  of  physiological  transport, 
molecular  behaviors  and  biochemical  reactions.  Parameter  estimation  mathematically  extracts  useful 
physiological  information  intrinsic  to  molecular  features  of  the  time  activity  curve  (TAC)  of  known 
normal  tissues  and  primary  malignancy  in  a  dynamic  PET  image  sequence.  These  procedures  have 
been  well-established  in  the  literature.  In  our  technique,  however,  once  the  parameter  estimations 
were  obtained  for  an  individual  patient,  the  corresponding  molecular  feature  spaces,  spanned  by  a  few 
molecular  feature  vectors  specifically  associated  to  this  patient,  was  formed  for  both  malignant  and 
normal  tissues.  A  feature-guided  filter  was  then  designed  for  this  patient,  based  on  his/her  own  feature 
spaces  rather  than  on  pooled-training  features,  considering  molecular  behaviors  may  vary  from  patient 
to  patient.  The  single  patient  tuned  filter  was  applied  to  the  dynamic  images  in  order  to  objectively 
suppress  the  background  tissue  interference  and  statistical  noise  from  the  contaminated  TACs  of 
small,  faint  metastatic  lesions  that  could  not  be  discerned  on  routine  visual  inspection.  The  filtered 
sets  under  two  hypothesized  data  models  were  processed  further  by  a  statistical  hypothesis  testing 
procedure.  This  procedure  made  hypotheses  on  the  data  according  to  tissue  types  (normal  or 
malignant)  and  then  tested  to  see  which  hypothesis  the  measured  data  were  most  likely  consistent 
with.  The  outputs  of  signal  detection  which  were  proportional  to  the  likelihood  (probability)  of  a 
given  hypothesis  being  true  or  false  for  each  measurement  were  used  as  an  adjunct  to  visual 
inspection. 


Figure  1 :  A  diagram  of  the  proposed  IDEMMF  system 
2.1.  Development  of  feature  extraction  and  analysis 

The  feature  extraction  and  analysis  in  dynamic  PET  images  have  been  studied  in  both  theory  and 
practice.  The  intensive  mathematical  analysis  has  been  performed  based  on  the  advanced  FDG-PET 
kinetic  modeling.  The  detailed  analysis  results  can  be  found  in  Appendix  A.  The  main  improvements 
on  the  existing  region  of  interest  based  techniques  are  briefly  reported  below. 
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2.1.1  TAC  feature  extraction 

Using  dynamic  phantom  data  with  known  ground  truth,  we  tested,  to  a  certain  degree,  how  the  time 
activity  curve  in  a  large,  visible  lesion  could  be  possibly  interfered  by  its  surrounding  background 
activity  (partial  volume  effect)  [1,2,  4,6]  through  the  currently  used,  imperfect  reconstruction 
techniques  and  how  the  accuracy  of  lesion  feature  related  parameters  might  be  affected  by  the 
common  ROI-averaged  time  activity  curve. 

We  have  performed  an  experimental  study  with  a  realistic  liver  phantom.  In  the  liver  phantom  three 
artificial  spherical  lesions  of  different  sizes  and  contrasts  were  placed  inserted  into  the.  Two  tracers,  nC 
and  !8F,  were  filled  into  the  liver  and  lesions,  respectively,  with  uniform  activity  distributions.  Dynamic 
data  were  acquired  of  the  phantom  with  ECAT953  2-D  whole  body  scanner.  This  allowed  the  radioactive 
decay  signature  of  the  two  radiopharmaceuticals  to  be  measured  as  pseudo-washout  data.  In  the  filtered 
backprojection  reconstructed  images,  the  smallest  lesion  with  7mm  interior  diameters  was  invisible  in  all 
dynamic  frames,  while  the  two  largest  lesions  can  be  clearly  visualized  in  the  FBP  images,  see  Figure  2 
(a)  -  (d)  for  an  illustration. 


fv 


••  I; 


«! 


\ 


% 


(a)  (b)  (c)  (d) 

Figure  2:  (a)  -(c)  The  last  frame  FBP  images  of  the  smallest,  median,  and  largest  lesions  in  the  liver  phantom;  (d)  The  time  activity 
curves  in  semi-log  scale:  true  mono-exponential  {unction  of  1SF  (red)  and  ROI  averaged  observation. 


The  largest  lesion  was  used  to  mimic  a  primary  tumor  detected  in  a  patient.  A  ROI  was  placed  in  it 
indicated  in  red  in  Figure  2  (c).  The  time  activity  curve  in  the  lesion  was  estimated  by  fitting  an  ROI 
averaged  observation.  The  resulting  curve  is  plotted  in  semi-log  scale  and  shown  in  Figure  2  (d).  As 
we  can  see  that  the  estimated  curve  in  blue  is  far  from  the  true  8F  time  activity  curve.  The  true  curve 
should  be  a  mono-exponential  function  and  a  straight  line  in  semi-log  scale.  Evidently,  the  estimated 
features  in  lesion  were  severely  contaminated  by  the  activities  of  nC  in  the  liver  background.  These 
features  can  not  adequately  characterize  lesions.  If  we  use  these  inaccurate  lesion  features  to  guide  a 
filter  design,  then  at  an  attempt  to  protect  the  lesion  features  during  the  filtering,  we  also  vulnerably 
keep  the  unwanted  background  interference.  This  will  definitely  affect  the  performance  of  filtering 
and  lower  the  SNR  gain  in  the  filtered  images. 

2.1.2  Non-invasive  blood  input  function  extraction 

Blood  input  function  is  required  in  molecular  feature  extraction  with  FDG  PET  time  activity  models 
[1,  8,9],  We  have  also  assessed  the  feasibility  of  replacing  invasive  blood  sampling  with  non-invasive 
blood  time  activity  curve  extracted  from  dynamic  image  data.  The  most  common  way  to  get  blood 
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function  is  invasive  blood  drawing;  however,  operationally  it  is  complicated,  time  consuming,  and 
difficult  for  patients  to  have  arterial  blood  draws  throughout  the  procedure.  A  critical  issue  with  the 
non-invasive  approaches  is  their  “ACCURACY”.  This  is  because,  besides  noise  interference,  the 
blood  function  is  often  interfered  also  by  surrounding  tissues.  There  are  two  non-invasive  approaches 
available:  one  is  a  simple  ROI  based  analysis  and  the  other  is  the  more  complex  FADS,  or  factor 
analysis  of  dynamic  structure  method.  FADS  is  well  known  for  its  capability  to  separate  time  activity 
curves  from  the  mixed  observations. 

Ten  lung  cancer  and  inflammation  patient  data  were  collected  for  the  study.  The  blood  samples  were 
invasively  collected  from  all  these  patients.  These  blood  samples  served  as  gold  standard  in  the 
evaluation. 


Figure  3:  (a)  The  last  frame  FBP  image  of  a  lung  cancer  patient;  (b)  The  ROI  used  to  non-invasively  extract  arterial  blood 
function;  (c)  The  three  blood  functions:  invasive  blood  samples  (red),  FADS  estimation  (green)  and  ROI-average  (blue). 


This  study  over  10  patients  concludes  that  when  compared  to  invasively  collected  blood  samples,  the 
blood  functions  estimated  from  FDAS  have  the  mean  square  errors  about  7  times  lower  than  those  of 
ROI  average  method.  FADS  approach  is  potential  to  outperform  ROI-based  method  for  replacing  the 
invasive  blood  sampling.  One  example  is  shown  in  Figure  3. 

FADS  processing  is  widely  used  in  cleaning  up  the  surrounding  background  time  activity 
interferences  in  the  blood  function  non-invasively  extracted  in  dynamic  images  at  an  area  with  strong 
arterial  activity.  Our  study  found  out  that  FADS  processed  blood  input  functions  has  much  more 
fidelity  to  the  invasively  collected  blood  samples  in  patients,  but,  on  the  other  hand,  we  have  also 
found  that  FADS  processing  failed  to  remove  the  interferences  superimposed  on  the  invisible  artificial 
lesions  in  both  experimental  and  digital  phantom  studies  due  to  low  lesion-to-background  ratios. 
Consider  the  detected  primary  tumor  and  the  normal  tissues  in  properly  selected  areas  are  usually 
strong  that  are  analogous  to  the  scenario  of  blood  function  extraction.  Thus,  we  will  resort  to  FADS 
processing  to  separate  the  unwanted  time  activities  contributed  from  surrounding  background  from 
the  observed  time  activities  in  the  primary  tumor  and  normal  tissues.  The  feature  parameters  will  be 
estimated  from  the  FADS  resulted  TACs,  instead  from  that  simply  averaged  in  ROI. 

2.2  Development  of  filtering  and  detection  algorithms 

Two  types  of  filtering  and  detection  algorithms  have  been  developed  which  exploit  the  differences 
between  normal  and  malignant  tissues  explored  in  the  feature  analysis.  The  details  are  reported  in 
Appendices  A  and  B  and  [8]  [9]. 
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2.3.  Evaluation  of  spatial-temporal  image  processing 
2.3.1  Evaluation  via  experimental  phantom 

Can  adding  temporal  information  to  spatial  processing  improve  SNR  in  small,  faint  or  invisible 
lesions?  This  is  a  critical  question  that  this  study  attempted  to  answer.  We  conducted  the  studies  with 
both  the  digital  (computer  generated)  phantom  and  the  liver  experimental  phantom  data  that  we 
believe  would  shed  light  upon  answering  the  question. 

The  above-mentioned  liver  phantom  shown  in  Figure  2  was  used  to  test  in  gaining  SNR  in  the 
smallest  lesion  through  spatial-temporal  processing.  Figure  4  (a)  -  (c)  shows  the  FBP  reconstructed 
images  of  this  lesion  at  Frame  1,15  and  23  and  Figure  4  (d)  is  the  MAP  reconstruction  of  Frame  23. 
This  confirms  that  the  lesion  is  invisible  in  all  the  frames  using  either  the  simple  or  advanced 
reconstruction  methods.  In  other  words,  by  exploiting  only  the  spatial  information  in  the  data  can  not 
improve  SNR  high  enough  to  make  the  lesion  become  detectable.  To  further  increase  SNR  in  the 
smallest  lesion,  making  a  use  of  the  intrinsic,  temporal  information  available  in  PET  dynamic  data  is 
pursuable.  We  first  tried  to  apply  FADS  to  separate  or  remove  the  liver  background  activities  and 
statistical  noise  superimposed  onto  the  lesion,  since  we  knew  that  FADS  is  well  credited  in 
decomposing  the  mixed  observation  from  visible  structures. 


Via  FADS  we  projected  the  dynamic  images  into  the  principle  components  of  the  images,  but  in  none 
of  the  decomposed  images  the  lesion  was  distinguishable.  In  Figure  4  (e),  the  first  two  principle 
component  images  are  presented  as  an  illustration.  This  indicates  that  a  direct  application  of  FADS 
provides  no  answer  to  early  metastasis  detection,  this  is  because  when  lesion  components  are  too 
weak  to  be  dominant  components,  the  principle  component  based  FADS  would  fail  to  separate  them 
from  the  noise.  Figure  4  (f)  is  the  output  result  of  the  metastasis  detection  shown  in  Figure  1.  The 
technical  details  are  given  in  Appendix  A.  In  this  test  data,  the  primary  tumor  features  were  extracted 
in  the  largest  artificial  lesion  shown  in  Figure  2  (c)  by  a  simple  ROI  average.  Although  the  accuracy 
of  ROI  based  feature  extraction  had  been  shown  not  high  enough,  the  results  of  this  study  reveal  that 
adding  temporal  information  to  image  processing  indeed  has  potential  to  increase  SNR  in  small 
lesions,  as  long  as  the  features  in  small  and  primary  lesions  are  similar  in  certain  degree. 


(e)  (f) 

Figure  4:  (a)-(c)  FBP  reconstruction  of  Frame  1,15  and  23;  (b)  MAP  reconstruction  of  Frame  23;  (c)  FADS 
decomposed  images;  (d)  filtered  image  using  the  simplest  version  of  proposed  methods. 


2.3.2  Evaluation  via  ROC  study  with  digital  phantom  study 

The  image  of  the  digital  phantom  is  shown  in  Figure  5  (a).  There  are  five  lesions  in  the  ellipse.  The 
lesions  were  assigned  with  the  time  activity  curve  collected  in  a  proven  lung  metastatic  malignant  and 
the  rest  of  the  ellipse  was  assigned  with  the  time  activity  curve  in  lung  normal  tissues,  respectively. 
These  curves  are  shown  in  Figure  5  (c).  By  forward  projecting,  the  dynamic  phantom  sinogram  data 
were  generated  and  certain  amount  of  Poisson  noise  was  added  to  the  projections.  The  noise  energy 
level  was  controlled  to  make  the  lesions  invisible  in  all  the  frames  of  FBP  reconstructed  images.  The 
last  frame  of  FBP  reconstructed  image  with  the  highest  SNR  is  presented  in  Figure  5  (b).  None  of  the 
five  lesions  can  be  visualized  as  desired. 

Assuming  that  the  time  activity  curves  in  the  lesions  and  the  background  are  known  in  prior,  we  tested 
at  each  pixel  whether  the  time  activity  curve  observed  is  more  similar  to  the  lesion  TAC  or  the 
background  TAC  in  the  sense  of  mean  square  errors.  If  at  a  pixel,  the  TAC  observation  is  closer  to  the 
lesion  TAC,  then  a  lesion  is  detected  at  that  pixel.  Otherwise,  the  pixel  is  declared  to  be  a  normal 
tissue.  50  sets  of  dynamic  phantom  data  were  generated  with/without  adding  digital  lesions  for  a  ROC 
study.  The  resulted  ROC  is  given  in  Figure  5  (d).  Although  in  reality  the  feature  knowledge  we  know 
about  lesion  and  normal  tissues  is  not  exactly  the  same  as  those  in  faint  or  invisible  metastases,  the 
resulting  ROC  demonstrates  the  potential  to  improve  early  metastasis  detection  from  almost  zero  to 
10% -20%. 


(a)  (b)  (c)  (d) 

Figure  5:  (a)  Digital  phantom  image;  (b)  FBP  reconstructed  image  of  (a);  (c)  Time  activity  curves 
assigned  to  lesions  (red)  and  background  (blue);  (d)  The  measured  ROC  curve 


2.3.3  Evaluation  via  animal  study 

We  evaluated  the  proposed  methods  for  early  detection  in  vivo  using  animal  study  with  a  miniature 
PET  scanner  (microPET,  Concorde  Microsystems,  Knoxville,  TN).  Primary  tumors  mimicking  the 
molecular  mechanisms  of  human  cancer  were  implanted  in  mice  subcutaneously. 

The  MDA-MB-435  human  breast  carcinoma  cell  line  was  used  to  grow  primary  tumors  in  the  study. 
During  the  period  of  tumor  formation,  dynamic  FDG  microPET  images  of  the  mice  were  acquired  to 
monitor  the  tumor  growth  at  early  stage  on  a  multiple-day  base.  Three  mammary  fat  mice  were  used. 
In  order  to  assure  the  early  stage  to  be  captured,  different  number  of  cells  (106,  0.5xl06  and  O.lxlO6) 
were  injected  into  the  same  mouse  under  the  two  arms.  In  this  means,  the  rates  of  tumor  growth 
should  be  different. 


9 


Dynamic  PET  imaging  of  FDG  (200  uci)  was  performed  on  a  MicroPET  R4  system  (Concorde 
Microsystems,  Inc).  Six  days  after  inoculation,  thirty-five  dynamic  data  frames  were  acquired  for  1  hr 
after  intravenous  injection  -  6xlsec,  4x3  sec,  10x30  sec,  5x60  sec,  and  10x300  sec.  Images  were 
reconstructed  with  the  OSEM  algorithm,  as  shown  in  Figure  6  (a).  The  corresponding  time-activity 
curves  of  the  two  tumors,  normal  tissues  and  heart  have  already  been  showed  in  Figure  6  (b).  The 
blood  input  function  was  measured  in  the  images  at  a  selected  ROI  in  the  heart  area.  Figure  6  (c)  is 
the  transformed  image  of  Figure  6  (a)  after  passing  through  the  constraint  temporal  filter. 
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Fig.  6.  Animal  study:  (a)  Reconstructed  FDG-PET  image  (6  days  after  inoculation);  (b)  TACs  drawn 
from  (a);  (c)  processed  image  of  (a)  by  using  the  constraint  temporal  filter;  (d)  Reconstructed  FDG- 
PET  image  (12  days  after  inoculation);  (e)  TACs  drawn  from  (d). 

We  also  present  another  reconstructed  image  in  Figure  6(d),  which  is  acquired  twelve  days  after 
inoculation.  Figure  6(e)  shows  the  corresponding  TACs  for  the  latterly  acquired  images.  Comparing 
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Figure  6(a)  and  Figure  6(b),  it  demonstrate  that  the  tumors  grow  bigger  and  the  slope  of  tumor’s  TAC 
becomes  steep  after  six  more  days.  Using  these  facts,  we  confirmed  that  the  findings  in  the  filtered 
output  image  shown  in  Figure  6(c)  were  correct. 

2.3.4  Evaluation  via  patient  data 

Patient  study 

Case  of  breast  cancer:  Fig.  7  (a)  shows  an  FBP  image  of  a  breast  cancer  data  archive  without 
attenuation  correction,  in  which  only  one  small  axillary  metastasis  can  be  visualized,  indicated  by  an 
arrow.  Attenuation  corrected  version  of  the  same  slice  is  shown  in  Fig.  7(b).  Applying  the  primary 
tumor  guided  filtering,  a  method  developed  during  this  project,  to  the  attenuation  corrected  dynamic 
images  identified  one  extra  metastasis  as  shown  in  Fig.  7(c).  These  positive  findings  have  been 
confirmed  by  a  follow-up  study  acquired  1.5  years  later  (the  patient  had  no  interval  treatment  of  her 
primary  or  metastatic  disease).  See  Figure  7  (d).  Between  the  two  studies  no  treatment  was  given  and 
the  disease  progressed.  Applying  a  threshold  to  the  filtered  image  in  Figure  7(c),  the  resulted 
detections  are  shown  in  Figure  7(e).  Overlaying  the  filtered  image  in  red  to  the  follow-up  image, 
Figure  7(f)  demonstrates  the  computer  identified  metastases  well  aligned  with  those  in  the  follow-up 
image.  The  slight  misalignment  between  lesions  and  detections  could  be  mitigated  if  an  advancing 
registration  accounting  for  non-rigid  transformation  was  applied.  Note  that  besides  the  two  confirmed 
lesions,  there  are  also  some  other  pixels  passing  the  threshold.  This  is  because  we  only  used  a  simple 
single  pixel  decision  criterion.  However,  if  we  add  the  spatial  size  of  lesions  to  be  detected  into  the 
decision  criterion,  then  we  can  get  rid  of  the  center  parts  in  Figure  7(e).  Moreover,  the  computer¬ 
generated  likelihood  of  malignance  will  be  reviewed  by  visual  inspection.  Combined  with  the 
additional  knowledge  of  the  patient’s  disease,  the  false  positive  findings  in  Figure  7(e)  can  be 
eliminated. 


(d)  (e)  (f) 

Figure  7:  (a)  FBP  image  of  a  breast  cancer  without  attenuation  correction;  (b)  Attenuation  correction  of  (a); 
(c)  the  filtered  image  of  (b);  (d)  the  follow-up  image  after  1.5  years;  (e)  the  detection  by  the  computer 
observer;  (f)  the  image  of  overlaying  (e)  on  (d). 
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Cases  of  lung  cancer  and  inflammation:  The  developed  spatial  temporal  filtering  guided  with  the 
features  extracted  on  a  patient  by  patient  base  has  also  been  applied  to  differentiate  between 
malignancy  and  inflammation.  Figures  8(a)  and  (b)  show  a  case  of  lung  cancer  and  the  corresponding 
filtering  output  image.  Evidently  the  lesion  was  enhanced  dramatically.  As  a  contrast,  Figure  9 
presents  a  patient  with  lung  cancer  and  joint  inflammation  where  the  joint  is  hot  in  the  late  frames  of  a 
dynamic  PET  scan.  The  simple,  single  template  defined  filter  (mentioned  in  last  annual  report)  was 
applied  to  process  the  data.  We  designed  the  filter  using  the  physiological  features  extracted  in  a  lung 
malignancy  shown  in  Figure  8(a).  In  the  filtered  image,  the  “hot”  spot  caused  by  benign  joint 
inflammation  was  filtered  out  the  as  normal  tissue  features.  The  result  is  shown  in  Fig.  9(b).  This  case 
demonstrates  that  the  “hot”  inflammation  in  FDGPET  images  may  mislead  visual  inspection  of  static 
image,  but  with  the  feature  guided  spatial-temporal  filtering,  we  can  have  a  better  chance  to  correctly 
distinct  it  from  malignances. 


(a)  (b) 

Figure  8:  A  lung  lesion 


(a)  (b) 

Fig.  9:  (a)  Joint  inflammation  in  a  patient  with  lung  cancer,  and  (b)  Filtered 
dvnamic  image  of  the  same  slice 


Remarks: 

The  three  year  grant  was  originally  funded  until  July  31,  2002.  However,  after  the  first  year’s  annual 
report,  it  was  brought  to  our  attention  that  we  needed  to  obtain  approval  from  both  the  U.S.  Army 
Medical  Research  and  Material  Command  Institutional  Review  Board  (US  Army  IRB)  for  patient 
data  collection,  because  although  patients  whose  data  to  be  used  in  this  study  will  have  a  PET  scan 
regardless  of  their  participation  in  this  study,  they  still  must  be  prospectively  recruited  for 
“additional”  (temporal-based)  images.  Since  November  2000  the  Army  IRB  had  have  a  hold  on  the 
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acquisition  of  the  clinical  data  under  this  grant  till  May  2002;  therefore,  we  had  not  obtained  all  of  the 
data  needed  to  complete  the  funded  objectives. 


We  have  requested  a  24  month  no-cost  extension  of  the  funding  period  for  this  project.  As  the  US 
Army  Medical  Research  Command  was  well  aware  of  the  delays  with  the  clinical  data,  the  request 
was  proved.  But  due  to  a  new  PET/CT  scanner  purchased  at  the  USC  PET  Clinical  Center,  the  old 
PET  scanner  which  has  the  capability  to  acquire  dynamic  scans  had  to  be  moved  to  and  installed  at 
our  research  facility,  thus  the  task  of  acquiring  75  patients  for  performance  analysis  and  evaluation  for 
this  project  was  not  fully  accomplished.  For  compensation,  we  added  three  groups  of  animal  studies, 
each  with  5  mice. 

The  major  accomplishments  of  this  activity  are  presented  as  follows: 


3.  Key  Research  Accomplishments 

The  main  accomplishments  of  the  project  are 

1.  Explored  temporal  physiological  differences  in  malignant  and  normal  tissues  based  on  the 
advanced  FDG-PET  kinetic  modeling;  mathematically  mapped  the  physiological  differences 
in  temporal  domain  onto  the  kinetic  (macro)  parameter  domain;  revealed  and  analyzed  such 
temporal  or  parametric  differences  in  frequency  domain  and  time-frequency  domain; 

2.  Assessed  the  accuracy  of  ROI-based  molecular  feature  extraction  techniques  with  the  real  liver 
phantom  data;  developed  the  factor  analysis  aided  feature  extraction  method  to  improve  the 
accuracy  of  feature  parameter  estimation;  evaluated  the  non-invasive  blood  input  function 
extraction  with  clinically  collected  blood  samples  at  the  USC  PET  center; 

3.  Designed  two  space-temporal  filtering  and  detection  criteria  to  identify  the  early  metastases 
embedded  from  the  unwanted  noise  and  background  interference  which  exploit  the 
physiological  features  in  malignant  and  normal  tissue  observed  in  FDG-PET  images; 

4.  Integrated  the  accomplishments  and  findings  listed  in  Items  1-3  into  a  software  prototype  of 
Intelligent  Detection  of  Early  Metastasized  Molecular  Feature  (IDEMMF)  system,  shown  in 
Figure  1; 

5.  Assessed  the  IDEMMF  system  with  digital  or  numerical  phantom,  experimental  phantom, 
animal  study  and  clinical  study. 
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4.  List  of  Reportable  Outcomes: 

4.1  Publications 

1.  C.  C.  Huang  and  X.  Yu,  "A  New  Method  of  Computer-Aided  Feature  Identification  for 
Lesion  Detection  in  PET-FDG  Dynamic  Study",  IEEE  Medical  Imaging  Conference  ,  Seattle, 
WA,  November  1999. 

2.  X.  Yu,  C.  C.  Huang  and  P.  S.  Conti,  "Comparison  of  Kinetic  Features  Extracted  in  OSEM 
and  FBP  Reconstructed  Dynamic  PET  Images  for  Oncology",  SNM  Annual  Conference,  San 
Louise,  Mo.,  June  2000. 

3.  X.  Yu  and  C.  C.  Huang,  "A  Fast  Method  to  Compute  Covariance  Matrix  in  Positron  Emission 
Tomography  Images",  IEEE  Medical  Imaging  Conference  2000. 

4.  C.  C.  Huang,  "Computer-Aided  Lesion  Detection  in  Positron  Emission  Tomography:  A 
Signal  Subspace  Fitting  Approach",  Ph.  D.  Thesis,  submitted  to  Electrical  Engineering 
Department  of  USC,  May  2001. 

5.  X.  Yu,  C.  C.  Huang  and  P.  S.  Conti,  "Assessment  of  ROI-based  Time  Activity  Analyses  in 
Dynamic  PET  For  Oncology  "  SNM  Annual  Conference  2001,  Toronto,  Ca.,  June  2001. 


6.  X.  Yu,  Z.  Li,  H.  Jadvar  and  P.  S.  Conti,  "Assessment  of  Non-invasive  Blood  Time  Activity 
Extraction  in  Dynamic  PET  Oncology”,  SNM  Annual  Conference  2002,  Los  Angeles,  CA., 
June  2002. 

7.  X.  Yu,  C.  C.  Huang,  and  P.  S.  Conti,  "Computer-aided  Metastasis  Detection  with  Dynamic 
PET”,  presented  at  the  Era  of  Hope  Conference,  Orlando,  FL.,  September  2002. 


8.  X.  Yu,  Z.  Li,  H.  Jadvar  and  P.  S.  Conti,  “Identification  of  Malignant  and  Benign  Lesions  in 
Dynamic  PET  Oncology”,  SNM  Annual  Conference  2003,  New  Orleans,  San  Louise,  June 
2003. 

9.  H.  Jadvar,  JR  Bading  and  X  Yu,  PS  Conti,  “Dynamic  FDG  PET  Kinetic  Analysis  of 
Inflammation  and  Cancer:  Preliminary  Results”,  SNM  Annual  Conference  2003,  New 
Orleans,  San  Louise,  June  2003. 

10.  J.  Chen,  X.  Yu,  “Enhanced  Dynamic  FDG-PET  Tumor  Detection  with  Constrained  Temporal 
Filtering”,  Proceedings  of  IEEE  International  Conference  on  Medical  Imaging,  Portland,  Or. 
October,  2003. 

11.  Z.  Li,  X.  Yu,  “Exploring  Frequency  Differences  Of  Physiological  Processes  To  Enhance 
Dynamic  FDG-PET”,  Proceedings  of  IEEE  International  Conference  on  Medical  Imaging, 
Portland,  Or.  October,  2003. 
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12.  J.  Chen,  X.  Yu,  “Rapid  Assessment  of  PET  Dynamic  Images  Using  Computer  Observers”, 
Proceedings  of  IEEE  International  Symposium  on  Medical  Imaging,  October,  2004. 

13.  Z.  Li,  X.  Yu,  “Computer  Aided  Lesion  Detection  with  Multi-channel  Time-Frequency 
Analysis”,  Proceedings  of  IEEE  International  Conference  on  Medical  Imaging,  October, 
2004. 

14.  X.  Yu  and  I.  S.  Reed,  “Theory  and  Algorithms  of  Rank  Reduction  for  Subspace  Filtering”, 
amended  submission  to  IEEE  Trans.  On  Information  Theory,  December  2004 

15.  J.  Chen  and  X.  Yu,  “Space-temporal  Analysis  and  Processing  of  Dynamic  PET  Images”, 
submitted  to  IEEE  Trans.  On  Information  Technology  in  Biomedicine. 

16.  Z.  Li  and  X.  Yu,  “Exploring  Physiological  Differences  in  Time-Frequency  Domain  to 
Improve  Tumor  Detectability  for  Dynamic  PET”,  submitted  to  IEEE  Trans.  On  Information 
Technology  in  Biomedicine. 

17.  Z.  Li,  X.  Yu,  “Computer-aided  Early  Metastasis  Detection  in  Dynamic  Positron  Emission 
Tomography”,  in  preparation  for  submission  to  IEEE  Trans.  On  Medical  Imaging  Processing. 


4.2  Education  Training  and  Graduation 
Master  Degree: 

1.  Chandhrasri,  Swias,  Electrical  Engineering,  January  1998  -  May  2000 

2.  Tian,  Ding,  Electrical  Engineering,  September,  2000  -  May  2002 

3.  Zheng,  Jianchang,  Computer  Science,  January  1996  -  December  1997 

Ph.  D.  Degree: 

1.  Huang,  ChungChie,  Electrical  Engineering,  September  1995  -  December  1999 

2.  Thanyasrisung,  Piyapong,  Electrical  Engineering,  May  1996  -  December  1997 

3.  Hu,  ChiaChang,  Electrical  Engineering,  September  1996  -  May  1998 

4.  Saghari,  Poorya,  Electrical  Engineering,  September  2000  -  May  2002 

5.  Chen,  Jiansong,  Electrical  Engineering,  September  2000  -  Present 

6.  Li,  Zheng,  Electrical  Engineering,  January  2001  -  Present 


5.  Conclusion 

The  goal  of  this  project  is  to  improve  detection  of  metastatic  axillary  breast  cancer  through 
sophisticated  physiological  modeling  and  statistical  signal  processing  techniques.  The  major  focus  of 
this  project  was  to  explore  temporal  physiological  differences  in  malignant  and  normal  tissues  based 
on  the  advanced  FDG-PET  kinetic  modeling  assessment;  assess  and  improve  the  accuracy  of  ROI- 
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based  molecular  feature  extraction  techniques  from  the  known  primary  breast  tumor;  design  the 
space-temporal  filtering  and  detection  criteria  to  identify  the  early  metastases  from  sever  background 
interference  and  count  noise;  integrate  the  developed  feature  extraction  and  filtering/detection  criteria 
into  a  software  prototype  of  Intelligent  Detection  of  Early  Metastasized  Molecular  Feature 
(IDEMMF)  system;  and  test  and  evaluate  the  prototype  with  phantom,  animal  study  and  clinical 
patient  study.  Our  theoretical  findings  include  mathematically  map  of  the  physiological  differences  in 
temporal  domain  onto  the  kinetic  (macro)  parameter  domain;  revealing  and  characterization  of  the 
temporal  or  parametric  domain  differences  in  frequency  domain  and  time-frequency.  The  evaluations 
on  a  small  scale  of  animal  and  patient  data  show  that  the  IDEMME  system  can  significantly  enhance 
the  metastatic  lesion  detection  by  exploring  the  temporal  differences  in  dynamic  FDG-PET  images. 
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Abstract 


As  with  all  nuclear  medicine  imaging  technologies,  lesion  detection  with  positron  emission  tomogra¬ 
phy  (PET)  scanning  of  2- [fluorine- 1 8] -fluoro-2-deoxy-D-glucose  (FDG)-labeled  tissues  is  restricted 
by  a  relatively  limited  spatial  resolution  and  a  low  signal-to-noise  ratio,  both  rendering  diagnosis 
by  visual  inspection  difficult  and  potentially  inaccurate.  In  this  paper,  a  computer-aided  lesion  de¬ 
tection  algorithm  combined  with  the  kinetic  time  activity  curve  (TAC)  feature  identification  from 
the  FDG-PET  dynamic  images  is  developed  to  assist  visual  inspection  for  lesion  detection.  We 
propose  that  the  unique  temporal  and  frequency  TAC  features  available  in  dynamic  FDG-PET 
images  can  be  distinguishedly  identified  for  normal  and  malignant  tissues,  then  be  usefully  applied 
in  statistical  hypothesis  tests  to  improve  the  detection  of  small  lesions.  The  hypothesis  is  that, 
the  TACs  extracted  from  FDG-PET  dynamic  images  may  have  variant  shapes  among  tissues  of 
different  size  and  location  within  a  given  organ,  but  the  sets  of  base  functions  weighted  to  compose 
the  TACs  which  could  be  physiological  process  functions  or  wavelets  base  functions  are  invariant. 
Thus,  our  approach  is  based  on  using  the  linear  invariant  subspaces,  defined  as  a  span  of  a  set  of 
base  functions,  to  model  the  lesion  and  normal  tissue  TAC  features.  The  lesion  and  normal  tissue 
subspaces,  estimated  from  the  known  tissue  types  of  data  by  the  signal-subspace-fitting  methods, 
are  used  for  unknown  lesion  detection  by  applying  a  matched  subspace  detection  criterion  derived 
from  a  generalized  likelihood  ratio  test.  The  experimental  results  using  FDG-PET  dynamic  studies 
of  cancer  patients  demonstrate  that  the  identified  TAC  features  are  not  only  distinct  and  invariant 
for  different  tissue  types  of  data,  but  also  useful  for  discriminating  lesions  from  normal  tissues  by 
incorporating  a  statistical  hypothesis  test.  To  quantitatively  compare  the  sensitivity  and  specificity 
for  different  computer  observers,  a  Monte  Carol  computer  simulation  of  dynamic  phantom  images 
are  simulated  for  achieving  the  receiver  operating  characteristic  analysis. 


1  Introduction 


Positron  emission  tomography  (PET),  a  rapidly  developing  technique  for  producing  in  vivo  radio- 
tracer  distribution  images,  provides  important  physiological  information.  Non-invasive  diagnostic 
procedures  using  PET  imaging  not  only  enable  physicians  to  assess  the  metabolic  activity  of  lesions 
in  vivo,  but  also  allow  patients  to  avoid  multiple  costly  procedures  and  to  benefit  from  guided  sur¬ 
gical  intervention.  PET  scanning  with  2-[fluorine-18]-fluoro-2-deoxy-D-glucose  (FDG)  is  improving 
the  diagnosis,  staging,  and  treatment  monitoring  of  a  variety  of  human  tumors  [3,  22]  and  is  becom¬ 
ing  useful  for  the  non-invasive  diagnosis  of  cancer.  However,  as  with  all  nuclear  medicine  imaging 
technologies,  lesion  detection  with  PET-FDG  images  is  restricted  by  a  relatively  limited  spatial  res¬ 
olution  and  a  low  signal-to-noise  ratio,  both  rendering  diagnosis  by  visual  inspection  difficult  and 
potentially  inaccurate,  especially  when  the  lesion  diameter  is  small.  Hence,  computer-aided  detec¬ 
tion  (CAD)  algorithms,  e.g.,  feature  identification,  have  been  developed  to  assist  visual  inspection 
in  PET  lesion  detection.  The  difficulty  is  that  spatial  features  such  as  the  shape  or  contrast  of  a 
lesion  in  an  image  are  often  varied  and  difficult  to  identify  [7,  8,  25].  Therefore,  the  PET  dynamic 
data,  providing  the  physiological  information,  becomes  straightforward  for  feature  identification 
[9], 

The  PET  dynamic  studies  have  been  widely  applied  in  medical  images  to  estimate  quantities, 
such  as  blood  volume,  flow,  and  local  cerebral  metabolic  rates  of  glucose  (LCMRGlc),  etc  [21, 
16].  For  example,  by  tracking  how  much  glucose  is  metabolized  in  different  areas  of  the  body, 
PET-FDG  imaging  enables  physicians  to  map  glucose  utilization.  The  compartmental  models  are 
usually  used  to  quantitatively  describe  the  tracer  kinetic  behaviors.  Physiological  factor  analysis 
has  been  applied  to  the  PET  dynamic  images  to  identify  fundamental  functions  useful  in  the 
compartmental  modeling  of  PET  data  [5].  Factor  analysis  of  dynamic  structures  (FADS)  was 
used  to  estimate  the  parameters  of  a  compartment  model  [1].  Because  FADS  is  used  to  solve  an 
under-determined  problem,  there  are  an  infinite  number  of  sets  of  factors  for  possible  solutions.  A 
procedure  called  information-based  factor  analysis  in  dynamic  structures  (IBFADS)  [14,  15]  was 
introduced  to  incorporate  a  priori  physiological  information  to  reduce  the  error  in  the  estimation 
of  correct  model.  Information  was  added  by  using  suitable  mathematical  models  to  describe  the 
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underlying  physiological  processes.  Principal  component  analysis  (PCA)  [17]  was  introduced  to 
enhance  clinically  interesting  information  in  a  PET  dynamic  imaging  sequence  in  the  first  few 
principal  component  images,  but  it  is  a  data-driving  technique  which  can  not  separate  signals  from 
high  noise  level.  A  mixture  model  was  introduced  to  estimate  the  parametric  images  with  improved 
mean  square  error  performance  and  low  computational  complexity  [16].  The  purpose  of  this  paper  is 
to  extend  the  use  the  PET-FDG  dynamic  study  for  lesion  detection  which  is  based  on  a  statistical 
parametric  imaging.  We  propose  that  the  unique  spatial  and  temporal  metabolic  time  activity 
curve  (TAC)  features  available  from  PET-FDG  dynamic  images  can  not  only  be  distinguishedly 
identified  for  normal  and  malignant  tissues,  but  also  be  useful  for  improving  small  lesion  detection 
by  incorporating  a  statistical  hypothesis  test. 

Prom  a  physiological  compartment  model  analysis,  the  PET-FDG  dynamic  study  for  a  given 
tissue  type  can  be  interpreted  by  a  homogeneous  compartment  model  which  is  governed  by  a  set 
of  kinetic  parameters  describing  the  passage  of  administered  tracer  for  this  given  tissue  type.  A 
homogeneous  compartment  is  associated  with  a  particular  dynamic  structure  called  a  physiological 
factor,  a  fundamental  TAC,  or  a  subspace  basis  [16].  A  PET  image  pixel  is  usually  modeled  to 
contain  more  than  one  physiological  factor  and  is  called  a  mixed  pixel  or  a  heterogeneous  pixel. 
It  is  assumed  that  a  heterogeneous  pixel  TAC  is  a  mixture  of  a  finite  number  of  independent 
homogeneous  compartments,  hence,  a  pixel  TAC  is  composed  of  a  combination  of  the  physiological 
factors  which  constitute  a  subspace.  In  this  paper,  we  hypothesize  that  each  (heterogeneous) 
tissue’s  kinetic  behavior  is  governed  by  its  own  set  of  kinetic  parameters,  which  suggests  that 
TAC  shapes  may  vary  among  different  tumors  within  a  given  patient,  but  will  occupy  the  same 
subspace  provided  they  have  the  same  kinetic  model  parameters.  For  example,  malignancies  can  be 
distinguished  from  normal  tissues  on  the  basis  of  biologically  determined  radiotracer  accumulation 
or  loss  rate.  Because  cancer  cells  divides  rapidly,  they  metabolize  glucose  at  a  much  higher  rate 
than  do  most  normal  tissues,  hence,  tumors  uptake  FDG  more  avidly  and  their  TACs  often  increase 
with  time  [11].  Therefore,  it  may  be  possible  to  use  the  invariant  subspaces  identified  from  visible, 
large  lesions  to  confirm  suspected,  but  not  unequivocally  identifiable,  small  lesions. 

We  use  the  linear  invariant  subspaces  to  model  lesion  and  normal  tissue  TAC  features.  By 
forming  the  PET-FDG  dynamic  data  into  a  temporal-spatial  matrix,  the  lesion  and  normal  tissue 
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subspaces  can  be  simply  estimated  by  using  the  signal-subspace-fitting  methods.  The  least  squares 
error  (LSE)  algorithm  is  used  for  parametric  subspace  estimation;  while  the  principle  component 
decomposition  is  used  for  non-parametric  subspace  estimation.  In  order  to  utilize  the  identified 
subspaces,  several  statistical  hypothesis  tests  are  derived  for  unknown  lesion  detection.  Both  su¬ 
perimposed  and  replacement  data  models  are  used  to  characterize  the  observations  of  large  lesions 
and  small  lesions.  Based  on  these  two  hypothesized  data  model,  the  statistical  tests  are  called  the 
matched  subspace  detectors  [19,  24]  which  can  be  derived  from  the  generalized  likelihood  ratio  test 
(GLRT)  principle.  The  matched  subspace  detector  is  actually  the  extension  of  the  rank-1  matched 
filter  to  the  multi-rank  filters.  In  this  paper,  both  the  rank-1  and  the  multi-rank  GLRTs  are  de¬ 
rived  for  the  lesion  detection.  The  rank-1  GLRT  was  derived  in  [19]  by  assuming  that  pixels  are 
uncorrelated  spatially.  Considering  that  image  pixel  values  are  not  statistically  uncorrelated  be¬ 
cause  each  count  contributes  to  all  pixel  values  [2],  the  rank-1  GLRT  needs  to  add  a  pre- whitening 
procedure.  The  spatial  covariance  matrix  required  for  the  pre-whitening  can  be  computed  using 
the  methods  described  in  [2,  10].  In  multi-rank  GLRT  approach,  it  is  assumed  that  the  spatial 
inter-pixel  correlation  in  each  frame  has  the  same  structure  but  different  energy  level. 

In  this  paper,  we  will  concentrate  on  the  filtered  backProjection  (FBP)  reconstructed  images. 
The  results  from  clinical  PET-FDG  dynamic  studies  show  that  the  physiological  features  estimated 
from  known  tumors  and  normal  tissues  are  distinct  and  invariant,  and  are  valuable  for  improving 
lesion  diagnosis  by  incorporating  a  statistical  test.  Finally,  to  quantitatively  compare  the  sensi¬ 
tivity  and  sufficiency  of  different  estimations,  detections  as  well  as  reconstructions,  a  Monte  Carol 
computer  simulation  of  phantom  dynamic  images  are  simulated  for  achieving  the  receiver  operating 
characteristic  (ROC)  analysis. 

In  Section  II,  the  PET-FDG  dynamic  data  modeling  described  by  a  three-compartment  FDG 
model  is  used  to  define  the  FDG  TAC.  Then  the  subspace  formation  is  introduced  based  on  a 
heterogeneous  pixel  assumption.  Section  III  presents  the  signal-subspace-fitting  methods  for  the 
lesion  and  normal  tissue  subspace  identification.  A  subspace  refining  method  based  on  a  subspace 
angle  criterion  is  shown  for  improving  the  subspace  separation.  Section  IV  adapts  and  extends 
the  GLRT  for  different  assumptions  of  pixel  correlation.  The  results  of  TAC  feature  identifications 
and  statistical  tests  for  lesion  detection  in  both  clinic  and  phantom  dynamic  data  are  presented  in 
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Figure  1:  4-k  three  compartment  model.  The  left  compartment  is  considered  as  the  blood  pool 
region  and  represents  the  concentration  of  FDG  in  the  plasma.  The  middle  compartment  contains 
the  FDG  in  tissue.  The  right  compartment  contains  the  phosphorylate  FDG  (FDG-6-P)  in  tissue. 


Section  V.  Section  IV  is  a  conclusion. 

2  Tracer  Kinetics  Modeling  for  Dynamic  FDG-PET 

In  this  paper,  we  will  focus  on  FDG-PET  dynamic  studies  because  FDG  is  widely  used  as  the 
tracer  in  clinical  PET  for  lung  and  breast  lesion  detection.  A  4-k  three  compartment  model,  one  of 
most  common  type  of  model,  is  used  to  model  the  kinetic  behavior  of  a  time  activity  curve(TAC). 
Two  cases,  namely  homogeneous  tissue  case  and  heterogeneous  tissue  case,  will  be  discussed. 

2.1  Homogeneous  pixel:  A  4-k  Compartment  Model 

Under  the  homogeneous  tissue  assumption,  the  FDG  tracer  distribution  in  the  tissue  can  be  modeled 
by  three  compartments  4-k  model  (Figure  2.1).  The  rate  constants  are  defined  as 

K\  =  the  transport  of  FDG  in  plasma  to  tissue, 

62  =  the  transport  of  FDG  in  tissue  to  plasma, 

63  =  the  phosphorylation  of  FDG  to  FDG-6-P, 
k\  =  the  dephosphorylation  of  FDG-6-P  to  FDG. 

where  *  indicates  decay-corrected  FDG  tracer  quantities. 

It  has  been  shown  [21]  that  the  radioactivity  for  a  homogeneous  tissue  Cus  is: 

Cusit)  =  — [(fc3*  +  k*4-  Pi)e~Plt  +  (/?2  -  fc3*  -  kl)e-^L]  0  Cp(t) 

P2  ~  PI 

=  [Afie-**  +  M2e-fc*]  ®  Cp{t),  (1) 
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where 


1  r  _ r 

A  =  2  (fc2  +  ^3  +  A*)  -  \Afc2  +  fcs  +  A*)2  -  4*2*3  , 

A  =  |  (A  +  K  +  kl)  +  ^ ( k*2  +  k*3  +  ki)*-  4 k*k*  , 

Mi  =  + 

Mi  =  (2) 

The  total  radioactivity  measured  by  the  PET  scanner  (under  noise- free  assumption),  (7(f), 
including  the  radioactivity  in  the  cerebral  blood,  can  be  represented  as 

C(t)  =  VbCp{t)  +  Ctis(t)  =  VbCp(t)  +  [. M ie-^f  +  M2e"^*]  ®  Cp(<),  (3) 

where  14  denotes  the  vascular  space  in  tissue  and  FDG  concentration  in  whole  blood  is  assumed 
to  be  VbCp{t)  at  time  t  [18]. 

2.2  A  Heterogeneous  pixel:  a  mixture  model 

Due  to  the  limited  spatial  resolution  of  the  PET  scanner  and  the  forward  and  backward  processes 
in  a  PET  system,  the  tissue  concentration  in  a  single  pixel  measured  by  the  scanner  is  the  mass- 
weighted  average  concentration  of  radioactivity  in  many  the  homogeneous  tissues  [12], [20].  Given 
that  a  mixed  (heterogeneous)  pixel  contains  J  homogeneous  subregions,  from  (3),  Cus(t),  the 
weighted  average  radioactivity  in  a  heterogeneous  tissue  can  be  written  as: 

J  J 

Ctis{t)  =  ^2  wiCtis  jit )  =  Y2  wj[Mije~Pljt  +  M2je~02it]  <8>  Cp(t),  (4) 

i=i  j= l 

where  wj  is  the  weighting  coefficient  for  the  jth  homogeneous  subregion.  Then,  the  total  radioac¬ 
tivity  measured  by  the  PET  scanner  in  a  heterogeneous  tissue,  (7(f),  can  be  represented  as: 

C(t)  =  VbCp(t)  +  Ctisit) 

3 

=  Vbcp(t)  +  J2  ^[Mye-A.i  +  M2je-^‘]  ®  Cp(t)  (5) 

j= 1 

Notice  that  Cp{t)  —  Cp(t )  ®  6(t )  ss  Cp(t )  <8>  when  (3  — >  oo.  For  real  data,  above  approxi¬ 
mation  works  well  when  (3  >  10.  By  this  way,  all  terms  in  (5)  are  expressed  as  convolution  of  blood 
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input  function  with  a  exponential  function.  Combining  some  constant  terms  and  rearranging  the 
subscript  j,  (5)  can  be  rewritten  as: 

J 

C(t)  =  ^2  aije-W  <g>  Cp(t)  (6) 

i= i 

where  atj  and  J  are  determined  by  M\j,  M2],  Wj,  14  and  J.  (6)  is  the  general  format  to  model  the 
TAC  measured  by  PET  scanner,  while  (3)  can  be  regarded  as  a  special  case  ( J  =  1)  for  it. 

In  all  models  listed  above  ,  the  plasma  input  function  Cp(t )  is  considered  as  s  given  function. 
There  are  several  techniques  available  to  obtain  this  function:  1)  normalized  population-based  input 
function,  2)  invasive  arterial  line  sampling,  3)  ROIs  of  aorta  or  large  vessels. 

3  Time-Frequency  Properties  of  TAC  and  Signal  Subspaces  for 
TAC 

The  FDG  malignant  (tumor)  cell,  in  general,  demonstrates  steady  increase  in  average  standardized 
uptake  values  (SUV),  while  SUV  in  normal  tissue  decrease  [29].  As  a  result,  the  observed  TAC  of 
malignant  tissue  increases  with  time,  whereas,  normal  tissue  TAC  washouts  with  time  [26].  This 
pattern  is  shown  to  be  true  in  at  least  80%  of  untreated  primary  cancers  [29].  Figure  4.3  illustrates 
this  temporal  difference  between  clinical  TACs  drawn  from  the  tumor  and  normal  ROI.  From  this 
starting  point,  followed  by  some  derivation  of  time-frequency  properties,  finally,  it  results  in  that 
the  tumor  and  normal  tissue’s  TAC  belongs  to  two  different  signal  subspaces. 

In  this  section,  we  will  first  exploit  the  temporal  difference  in  frequency  domain  to  characterize 
the  behaviors  of  the  individual  physiological  process  functions.  Then,  by  decomposing  the  TACs 
into  summation  of  physiological  process  functions  or  alternatively  decomposing  the  TACs  into 
summation  of  wavelet  base  functions,  the  time-frequency  differences  between  the  tumor  and  normal 
TACs  are  revealed.  This  differences  imply  that  tumor  and  normal  TACs  can  be  approximated  by 
different  base  function  sets,  which  means  the  tumor  and  normal  TACs  belong  to  different  signal 
subspaces. 
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3.1  Time  Domain  Properties  of  TAC 


Define  the  physiological  process  function  [27]  as: 

f(P,t)  =  Cp(t)®e-0t,t>O  (7) 

where  Cp(t)  is  the  arterial  blood  input  function,  (3  is  the  decay  coefficient.  Since  Cp(t)  is  same  for 
an  individual  patient,  it’s  f3  value  that  affect  the  temporal  properties  of  different  f(/3,t).  Before 
the  further  derivation,  two  constrains  are  added  to  Cp(t)  under  the  noise  free  assumption.  First, 
Cp(t)  must  be  non-negative,  i.e.  Cp(t)  >  0,t  >  0;  second,  there  exist  a  To  for  each  individual 
Cp(t)  such  that  dCp(t)/dt  >  0,  t  <  To  while  dCp(t)/dt  <0,  t  <  Tq  These  constrains  have  been 
verified  by  a  large  patient  popularity  [28]  as  well  as  by  plasma  TAC  model  [4].  With  these  two 
constrains,  we  proved  (Appendix)  that  for  a  given  time  T,  T  >  To,  there  exists  a  threshold  (3*  in 
dynamic  FDG-PET  for  each  individual  patient,  such  that  the  df(f3,  t) / dt\p-p* tt=T  —  0-  Then  for 
/ 3  value  which  satisfies  0  </?</?*,  df(/3,  t) / dt'r0<t<T  >  0,  while  for  [3  >  (3*,  df((3,  t)/dt\t^T  <  0. 
Remind  that  physiological  process  function  f{(3 ,  t)  for  tumor  cell  will  go  up  along  with  time  while 
the  normal  cell’s  physiological  process  function  f(/3,  t )  will  go  down  with  time.  Combining  with 
the  temporal  properties  we  just  drawn  above,  we  can  conclude  that  /(/?,  t)  for  tumor  cell  will  have 
(3  <  /?*;  ,  while  f{(3,  t )  for  normal  cell  will  have  (3  >  (3* . 

Using  (6)  and  (7),  the  TAC  measured  by  PET  scanner  in  the  heterogeneous  tissue  C(t),  can  be 
modeled  as: 

J  J 

C(t)  =  £  aje-**  ®  Cp(t)  =  Y,  t>  0  (8) 

j=i 

When  decomposing  a  TAC  into  the  linear  combination  of  many  physiological  process  functions, 
we  hypothesis  that  the  physiological  process  functions  with  (3  smaller  that  (3*  are  associated  to 
malignancy  tissue,  while  the  physiological  functions  with  [3  bigger  than  (3 *  are  related  to  normal 
tissue.  So,  the  physiological  process  functions  can  be  divided  into  two  sets:  {/(/?,  t)\{3  <  (3*}  and 
{/(/?,  t) |/3  >  /?*}.  The  first  set  forms  a  linear  subspace  in  which  tumor’s  TAC  lies,  the  second  set 
form  a  linear  subspace  in  which  normal  tissue’s  TAC  lies. 
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Figure  2:  A  bunch  of  the  physiological  process  functions  /(/?,  t)  when  (3  varies  from  0  to  1 


3.2  Frequency  Domain  Properties  of  TAC 

From  the  linear  system  point  of  view,  a  physiological  function  is  the  output  of  the  arterial  input 
function  Cp(t )  passing  through  a  linear  time  invariant  system  with  impulse  response: 

h((3,  t)  =  e_/3i,  f>0  (9) 

Applying  the  Fourier  Transform  to  the  both  sides  of  this  equation  (8): 

J 

C(Lj)  =  '52ajcP{u,)H(pj,u>)  (10) 

i=i 

where  Cp{lo)  is  the  Fourier  Transform  of  blood  function  Cp(t),  H(f3j,ui)  is  the  Fourier  Transform 
of  h(/3j,t ): 

H(J3,u)  =  l/(P  +  jcj),  \H((3,u)\  =  l/Vf32  +  u2  (11) 

Notice  that  the  smaller  the  (3,  the  lower  the  filter’s  bandwidth.  So,  for  a  certain  patient  (with 
a  given  blood  input  function),  the  physiological  process  function  with  smaller  (3  contains  less  high 
frequency  components  than  those  with  the  larger  (3.  Recall  our  hypothesis  described  in  the  previous 
subsection,  (10)  (1 1)  indicates  that  the  TAC  in  malignance  has  fewer  high  frequency  components 
than  that  of  normal  tissue. 
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Figure  3:  Two-stage  wavelets  analysis  filter  bank.  a;[n]  is  the  input  signal,  vt [n]  is  the  output 
wavelet  coefficients  of  channel  i.  h,L  represents  the  lowpass  subband  filter,  ha  represents  the 
highpass  subband  filter. 

3.3  Time-Frequency  Domain  Differences  between  Tumor  and  Normal  TAC 

In  this  section,  the  wavelet  transformation  will  be  employed  so  as  we  can  get  frequency  domain 
information  as  well  as  short  period  time  information.  Further  more,  the  wavelet  transformation 
also  provides  an  alternative  way  to  express  and  identify  the  tumor  and  normal  tissue’s  subspaces. 
The  discrete-time  wavelet  analysis  and  synthesis,  in  vector  format,  can  be  written  as: 


v  =  Tay 

(12) 

y  =  Tsv 

(13) 

where  y  is  the  input  signal,  v  is  the  wavelet  coefficients,  Ta  is  the  analysis  matrix  and  Ts  is  the 
synthesis  matrix  whose  columns  are  discrete-time  wavelet  base  functions: 

Ts  =  [ipi,^2,  ■■■  ,i>k]  (14) 

where  ipi,  i  =  l...k  is  the  discrete-time  wavelet  base  functions,  k  is  the  total  number  of  base  functions. 

In  practice,  the  discrete-time  wavelet  analysis  can  be  implemented  by  subband  linear  filter 
bank[31].  Suppose  hi, [n]  represents  the  low-pass  subband  filter,  hjj [n]  represents  the  high-pass 
subband  filter,  y[n]  represents  the  input  signal.  Figure  3.3  illustrates  a  two-stage  cascading  structure 
of  wavelet  analysis.  When  passing  a  TAC  through  a  multi-channel  filter  bank,  the  outputs  of 
different  channels  represent  the  components  in  different  frequency  bands.  The  higher  the  stage,  the 
lower  the  frequency.  The  outputs  at  the  same  channel  carry  the  time  information. 
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Figure  4:  Wavelet  tilings  of  noise  free  tumor  and  normal  tissue’s  TAC.  (a)  clinical  noisy  TAC  and 
denoised  TAC  for  normal  tissue.  The  denoise  is  done  by  model  base  least  square  curve  fitting 
method,  (b)  wavelet  tiling  of  denoised  normal  TAC.  The  X-axis  represents  the  time.  Y-axis 
represents  the  frequency  band  number,  i.e.  the  channel  number  in  the  figure  3.3.  The  higher  the 
channel  number  the  lower  the  frequency  band,  (c)  clinical  noisy  TAC  and  denoised  TAC  for  tumor, 
(d)  wavelet  tiling  of  denoised  tumor  TAC. 


Plotting  the  wavelet  coefficients  of  a  certain  signal  in  2D  format,  in  which  X  axis  represent 
the  time  and  Y  axis  represent  the  frequency,  we  will  get  time- frequency  ’’tiling”  which  shows 
the  time-frequency  distributions  of  that  signal.  Figure  3.3  illustrates  the  time-frequency  tiling  (the 
amplitudes  of  the  coefficients  are. shown)  of  clinical  normal  and  tumor  TACs  using  wavelet  analysis. 
It  is  shown  that  (i)  tumor  TAC  has  fewer  high  frequency  components  than  that  of  normal  tissue,  (ii) 
energy  of  tumor  TAC  increases  in  the  later  time  slots,  while  the  energy  of  normal  TAC  decreases. 
These  observed  patterns  of  time-frequency  tilings  concord  with  the  tumor  and  normal  TACs’  time 
domain  and  frequency  domain  properties  derived  before. 

The  purpose  of  involving  the  wavelet  analysis  here  is  to  express  the  subspaces  using  wavelet  base 
functions.  According  to  (13)  (14),  the  wavelet  transformation  can  be  regarded  as  the  decomposition 
of  the  input  signal  to  the  weighted  summation  of  independent  wavelet  base  functions.  The  examples 
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Figure  5:  The  Daubechies  4  wavelet  base  functions  for  512  points  seven  levels  filter  bank.  The  base 
functions  for  different  levels,  i.e.  for  different  frequency  bands,  are  shown. 

of  base  functions  are  shown  in  figure  3.3.  The  different  patterns  of  time-frequency  tilings  between 
the  noise  free  tumor  and  normal  TACs  imply  that  the  tumor  and  normal  TACs  can  be  approximately 
expanded  by  different  subset  of  wavelet  base  functions.  In  another  word,  the  tumor  and  normal 
TACs  belong  to  the  different  linear  subspaces. 

4  Subspace  Formation  and  Identification 

In  this  section,  we  first  give  a  general  formation  of  the  subspaces,  then  two  subspace  identification 
methods  are  provided.  The  clinical  studies  show  that  our  hypothesis  of  existence  of  tumor  and 
normal  tissue’s  subspaces  is  reasonable:  the  identified  subspaces  can  express  tumor  and  normal 
TACs  accurately.  Beside  the  two  sets  of  base  functions  will  be  discussed,  namely  physiology  pro¬ 
cess  functions  and  wavelet  base  functions,  the  cubic  B-splines,  which  are  used  in  4D  continuous 
reconstruction,  potentially  could  be  used  as  the  base  functions.  By  this  way,  it  prospectively  can 
link  the  reconstruction  process  with  tumor  detection  task. 

4.1  The  General  Formation  of  Subspaces 

The  TAC  received  by  the  PET  scanner  with  noise  can  be  modeled  as: 

OO 

y{t)  =  C(t)  +  n(t)  =  ^  afej(f)  +  n(t)  (15) 

i=l 
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where  {ei(t),  e2(f), ...}  is  a  certain  set  of  independent  continuous  base  functions  (without  loss  of 
generality,  infinity  number  of  base  functions  are  used  here.  For  real  data,  finite  number  of  base 
functions  should  be  accurate  enough).  After  sampling  this  continuous  signal  at  N  time  points,  it 
becomes: 

k 

y  =  aiei  +  n  =  Ea  +  n,  (16) 

i=  1 

where  y,  e*, n  €  RN,  e,,i  =  1, 2, ...,  k  are  independent;  a  =  [ai,  a?, ..., a*,]'  £  Rk;  E  =  [ei,  e2, ...,  e*.] 
is  a  N  x  k  matrix.  Because  {ei,e2, ...,  ej.}  is  a  set  of  base  function  defined  on  RN ,  k  <  N  should 
be  satisfied.  The  noise  n  is  assumed  to  be  additive  Gaussian  noise  in  this  paper. 

According  to  the  subspace  assumption  we  mentioned  before,  span{ei,e2,  ...,e/c}  is  the  union 
of  two  subspaces:  one  is  the  subspace  for  the  pure  tumor  TAC,  denoted  by  (H);  another  is  the 
subspace  for  the  pure  normal  TAC,  denoted  by  (S): 


(17) 


(H)  =  span{eil,ei2,...,eip},{i i,i2,....,ip}  C  {1,2,...,  A:} 

(S)  =  span{ejl,eh,...,ejq},{ji,j2,-;jq}  C  (1, 2, ...,  k} 
where  p  is  the  dimension  of  tumor’s  subspace  and  q  is  the  dimension  of  normal  tissue’s  subspace. 
(H)  and  (S)  may  have  some  base  functions  in  common,  but  each  of  the  subset  must  contain  some 
base  functions  which  are  not  included  in  another  subset.  Define  matrix  H  and  S  as: 


H  —  [  e . , ,  ej2 , . . . ,  ej 


S  —  [e jj ,  ej2 , . . . ,  e 


(18) 


■0q\ 


According  to  (16)  and  (18),  the  noisy  TAC  y  sampled  from  the  normal  tissue  region  and  tumor 
region  can  be  respectively  modeled  as: 


Ho  ( normal )  :  y  =  Sa  +  n 

Hi  (tumor)  :  y  =  Ha  +  n 

where  a  £  R q  for  normal  case,  and  a  £  itRP  for  tumor  case. 

4.2  Subspace  Identification  Methods 


(19) 

(20) 


Given  the  P  pixels  tumor  or  normal  tissue’s  region  of  interest  (ROI)  date,  the  subspace  identification 
problem  is  to  identify  H  or  S  respectively.  Using  (20)(20),  the  data  model  can  be  written  as,  for 
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P  pixels  tumor  tissue’s  ROI  data: 


Y  =  HA  +  N  (21) 

for  P  pixels  normal  tissue’s  ROI  data: 

Y  =  SA  +  N  (22) 

where  Y  =  [yi,y2)  — >yp],  A  —  [aiia2,  N  =  [ni,  112, n/>].  Because  the  identification 

problem  is  similar  for  tumor  and  normal  tissue,  only  the  case  for  tumor  case  is  discussed. 

4.2.1  Parametric  :  Least  Squares  Estimation  (LSE) 

For  the  tumor  ROI  data  described  in  (21),  the  subspace  fitting  using  the  LSE  is  defined  as  [23]: 

H,  A  =  argmin  ||Y  —  HA|||>,  (23) 

H,A 

where  ||  •  ||p  denotes  the  Frobenius  norm.  It  can  be  derived  [23]  that: 

H  =  arg  maxtr{PHYYTP£}  (24) 

H 

where  PH  =  H(HtH)-1Ht  is  the  orthogonal  projection  matrix  that  projects  onto  the  column 
space  of  H  and  “tr”  is  the  trace  operator.  When  the  measured  blood  function  Cp(t)  is  available, 

each  column  of  H  can  be  modeled  using  (6),  otherwise  (??)(6)  will  be  used.  By  this  way,  the 

only  unknowns  in  H  are  the  exponential  parameters  /3j,  so  the  LSE  becomes  to  search  the  fij  over 
a  reasonable  range,  i.e.  [0, 1]  in  this  case(One  f3j  is  corresponding  to  blood  function  and  so  as 
considered  a  known  term).  Newton-Raphson  [30]  or  spectral  analysis  method  [13]  can  be  used  to 
get  the  numerical  solution  for  this  problem. 

4.2.2  Non-Parametric:  Wavelet  Transformation  Method 

Wavelet  based  subspace  identification  procedure  includes  three  steps.  First,  using  synthesis  filter 
bank  (13)  to  compute  the  complete  set  of  discrete-time  wavelet  base  functions.  Second,  using 
analysis  filter  bank  (12)  to  compute  the  mean  value  of  wavelet  coefficients  y  for  TACs  at  the  given 
P  tumor  pixels.  Third,  the  base  functions  for  p  dimension  subspace  of  the  tumor  are  wavelet  base 
functions  with  the  first  p  most  significant  wavelet  coefficients.  So, 

H  1  V'^2  5  •  ■  • ,  ]  (25) 
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Figure  6:  Clinic  FBP  reconstructed  image  and  ROI  TACs.  (a)  the  last  frame(28tft)  of  the  PET- 
FDG  dynamic  FBP  reconstructed  images,  where  three  5x5  pixels  ROIs  are  chosen:  LI  in  the  large 
lesion  (primary),  L2  in  the  small  lesion  (metastasis)  and  BG  in  the  background  (normal),  (b)  the 
average  TACs  for  ROI  LI,  ROI  L2  and  ROI  BG  respectively. 

The  subscripts  *1, ...  ,  ip  are  chosen  such  that  vn , ...  ,Vjp  are  first  p  most  significant  values  in 
i  =  1  ...k. 

4.3  Subspace  Identification  Results 

A  clinical  28-frame  PET-FDG  dynamic  study  of  lung  cancer  is  shown  here  to  verify  our  hypothesis 
about  the  existence  of  the  subspaces  and  the  subspace  identification  methods  (The  detailed  infor¬ 
mation  about  the  clinical  data  is  listed  in  the  next  section).  There  are  two  lesions  in  this  study, 
clearly  visualized  large  tumor  (primary)  and  another  barely  seen  smaller  tumor  (metastasis).  Fig¬ 
ure  4.3  shows  the  clinic  image  and  ROI  TACs  sampled  from  the  image.  Our  clinical  study  shows 
the  identified  subspaces  are  accurate  and  separable. 

Parametric  LSE  Method:  The  normal  and  tumor  TAC’s  subspaces  estimated  by  the  LS  method 
were  found  to  be  spanned  by  (e_0lt,  e~02t,  e~8:it}  and  {1  —  e~/3lt,e~^2t,te~^st},  respectively,  and 
the  parameter  sets  were  {$i  =  0.00025,  62  =  0.005,  63  =  0.05}  and  {/?i  =  0.005,  /%  —  0.0223, 
/?3  =  0.00125}  for  BG  and  LI,  respectively.  Figure  (??)  show  the  curve  fitting  results  using  sub¬ 
spaces  identified  by  LSE.  Table  (4.3)  shows  the  curve  fitting  MSE.  Both  the  figure  and  MSE  show 
that  the  identified  tumor  subspace  H  can  fit  tumor  TAC  well  while  fit  normal  TAC  badly;  the 
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LARGE  TUMOR  LI  BACKGROUND  BG  SMALL  TUMOR  L2 


Figure  7:  Comparison  of  accuracy  and  separability  of  extracted  lesion  and  normal  tissue’s  subspaces 
using  least  square  estimation.  The  H  is  estimated  tumor’s  subspace  using  TAC  data  in  LI,  S  is 
estimated  normal  tissue’s  subspace  using  TAC  data  in  BG.  (left)  TAC  in  LI  represented  by  the 
subspace  S  and  TAC  in  LI  represented  by  the  subspace  H.  ( middle )  TAC  in  BG  represented  by 
the  subspace  H  and  TAC  in  BG  represented  by  the  subspace  S.  (right)  TAC  in  L2  represented  by 
the  subspace  H  and  TAC  in  L2  represented  by  the  subspace  S. 


identified  normal  subspace  S  can  fit  normal  TAC  well  while  fit  tumor  TAC  badly.  In  one  word,  the 
identified  subspaces  are  accurate  and  separable. 


Ll  TAC 

L2  TAC 

BG  TAC 

taae 

0.0197 

0.0123 

0.0226 

BEI 

0.0536 

0.0266 

0.0195 

Table  1:  Curving  fitting  MSE  for  different  TACs.  The  subspaces  <  H  >  and  <  S  >  are  estimated 
by  LSE 


Non- Parametric  Wavelet  Transformation  Method:  The  equations  (12)(??)  (25)  and  the  base 
functions  choosing  criteria  are  applied.  Setting  p  =  8,  q  =  12,  the  base  functions  for  tumor  and 
normal  subspaces  are  illustrated  in  Figure  4.3.  Figure  4.3  shows  the  curve  fitting  results  using 
subspaces  identified  by  wavelet  method.  Table  (4.3)  shows  the  curve  fitting  MSE. 

The  results  from  both  methods  show  that  the  estimated  subspaces  S  and  H  can  respectively 
represent  the  normal  and  lesion  TAC  very  accurately  (with  small  MSE).  At  the  same  time,  S  and 
H  are  separable:  S  can  not  represent  the  tumor  TAC  accurately  (with  large  MSE)  and  H  can  not 
represent  normal  TAC  accurately. 
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Subspace  <S>,  Normal 


Figure  8:  The  tumor  subspaces  (H)  and  normal  subspace (5')  extracted  using  wavelet  method 


LARGE  TUMOR  LI  BACKGROUND  BG  SMALL  TUMOR  L2 


Figure  9:  Comparison  of  accuracy  and  separability  of  extracted  lesion  and  normal  tissue’s  subspaces 
using  wavelet  method.  The  H  is  estimated  tumor’s  subspace  using  TAC  data  in  LI,  S  is  estimated 
normal  tissue’s  subspace  using  TAC  data  in  BG.  Notice  that  some  estimated  TAC  overlap  with  the 
patient  TAC.  (left)  TAC  in  LI  represented  by  the  subspace  S  and  TAC  in  LI  represented  by  the 
subspace  H.  ( middle )  TAC  in  BG  represented  by  the  subspace  H  and  TAC  in  BG  represented  by 
the  subspace  S.  (right)  TAC  in  L2  represented  by  the  subspace  S  and  TAC  in  L2  represented  by 
the  subspace  H. 


LI  TAC 

L2  TAC 

BG  TAC 

<  H  > 

0.0154 

0.0208 

0.0296 

<  S  > 

0.0887 

0.0446 

0.0134 

Table  2:  Curving  fitting  MSE  for  different  TACs.  The  subspaces  <  H  >  and  <  S  >  are  estimated 
by  wavelet  analysis  method 
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5  Matched  Subspace  Detection 


The  generalized  likelihood  ratio  test  (GLRT)  criterion  is  employed  to  solve  the  detection  problem. 
By  using  the  principle  of  the  GLRT,  the  matched  subspace  detectors  [19,  24]  will  be  derived  under 
two  different  assumptions  of  noise’s  time  correlation  and  pixel’s  spatial  correlation.  The  first  one 
assumes  that  the  pixel  is  uncorrelated  spatially  and  the  noise  is  white  in  time.  The  second  one 
assumes  the  pixel  is  spatially  correlated  and  noise  is  colored.  For  each  of  the  detector,  replacement 
data  model  and  superimposed  data  model,  axe  applied. 

5.1  Hypothesized  Data  Model 

Depending  on  modeling  whether  the  lesion  is  mixed  with  surrounding  normal  tissue  or  not,  two 
data  models  are  used  for  the  detection  hypothesis. 

Replacement  Model:  The  lesion  absence  hypothesis  Ho  says  that  the  data  consists  of  a  sum  of 
normal  tissue  signal  x0  and  the  noise  no-  The  lesion  existence  hypothesis  H\  says  that  the  data 
consist  of  a  sum  of  lesion  signal  xi  and  the  noise  ni.  That  is, 

Rb:y  =  xo  +  n0  and  #1  :y  =  xi+ni.  (26) 

where  xi  =  H0,  9  €  Hp\  xo  =  S<f>,  <f>  €  Hq.  no  and  ni  are  additive  Gaussian  noise. 

Superimposed  Model:  The  lesion  absence  hypothesis  Ho  is  same  as  replacement  model:  the  data 
consists  of  a  sum  of  normal  tissue  signal  xo  and  the  noise  no-  The  lesion  existence  hypothesis 
Hi  says  that  the  data  consist  of  a  sum  of  lesion  signal  xi ,  normal  tissue  signal  xoi  (which  means 
x0i  €  S,  but  xoi  may  not  be  equal  to  xo)  and  the  noise  ni.  That  is, 

H0:  y  =  x00  +  n0  and  Hi  :  y  =  (xi  +  x0i)  +  ni ,  (27) 

where  xi  =  H 6,  9  €  7 Zp;  xoo  =  S(f)o,  4>o  G  Hg;  xoi  =  S<f>i,  cj) \  G  lZg.  no  and  ni  are  additive 
Gaussian  noise. 

5.2  Generalized  Likelihood  Ratio  Test 

Given  hypothesis  test  with  model: 

H0:  y  =  x0  +  n0  and  #i:y  =  Xi+ni.  (28) 
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The  generalized  likelihood  ratio  test  is: 


Ky)  = 


*(*i|y) 

o|y) 


*8° 


(29) 


where  l(-)  =  ln(p(-)),  p(-)  is  the  probability  density  function,  xi  and  xo  are  maximum  likelihood 
estimate(MLE)  of  xi  and  xq  respectively. 


5.3  Matched  Subspace  Detectors 


Case  1:  White  Noise  Without  Spacial  Correlation 

For  the  replacement  model  defined  in  (26)  with  the  assumption  that  the  noise  nj  is  normal  with 
zero  mean  and  covariance  matrix  of  I  and  the  signal  x,  obey  the  linear  subspace  model  xj  =  HOj 
and  xq  =  S6>o-  Using  (29) [19]: 


h(y)  = 


J(Mily) 


A)-N/2 

V 


cp{- 


ni|l2+7^ll“o||2} 


oly)  #0  26>f 

Substituting  the  MLE  of  <jj,  6  and  0  into  it,  finally  result  in: 


2^o 


P(0,v  i) 

o) 


(30) 


£lM_  ll^olli  _  yTp5y 
(y) "  iifiiiil  "  --TT‘±- 


2  yTPi^y  </f°  m 


(31) 


where  P#  is  the  orthogonal  projection  onto  the  complement  subspace  of  (H)  and  Pjr  is  the  or¬ 
thogonal  projection  onto  the  complement  subspace  of  (S),  which  are  defined  as: 

Pg  —  I  -  S(STS)_1ST 
P#  =  I  —  H(HtH)"1Ht 

Similarly,  for  the  superimposed  model  defined  in  (27),  the  GLRT  becomes  [19]: 


ToM  =  H^oill  _  yTpsy  >hx 
2(y)  llfiilli  yTHsy<Ho 


(32) 


(33) 


where  P^s  is  the  orthogonal  projection  onto  the  complement  subspace  of  (HS)  =  (H)  (J(S),  which 
is  defined  as: 


Pjjs  =  I  -  HS(HStHS)-1HSt  (34) 

Case  2:  Colored  Noise  With  Spacial  Correlation 

It  is  assumed  that  the  spatial  inter-pixel  correlation  in  each  frame  has  the  same  structure  but 
different  energy  level.  Let  the  P-dimensional  vector  f,  =  [/,(!),  /,  (2),  •  •  •  ,  fi(P)]T  denote  the  spatial 
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pixel  data  for  the  z-th  frame  and  denote  F  =  [fi,  (2,  •  •  ■  ,  fjv].  We  assume  that  all  the  columns  of 
the  matrix  F  are  independent.  Hence  Y  =  [yi,  y2,  •  •  •  ,  y p]  =  FT.  The  derivation  of  the  multi-rank 
GLRT  for  the  constrained  correlation  structure  is  based  on  a  replacement  model,  i.e., 


r  Ho  :Y  =  Xq  +  N0, 
\  Hi  :  Y  =  X1+N1. 


The  GLRT  can  be  derived  as  : 


1  +  (G~1/2Rs)TPg  (G-^Rs) 

31  ’  l  +  (G-1/2Rs)Tp^(G-l/2R§)  <H°  m 


(36) 


where  R  is  the  covariance  matrix,  s  =  M  1/21,  s  =  (s~s)1//2s  ,G  =  RRT  —  (Rs)(Rs)T.  The 
matrix  M  is  the  correlation  structure  for  each  frame. 


6  Result  and  Discussion 

Three  categories  of  dynamic  PET  studies  are  involved  to  demonstrate  our  methods.  Physical 
phantom  study  was  designed  to  validate  the  efficacy  of  MSD  for  dynamic  PET  lesion  detection. 
Digital  phantom  study  and  clinical  study  were  designed  to  validate  the  subspace  identification 
methods  as  well  as  MSD  with  estimated  subspaces. 

6.1  Dynamic  PET  Physical  Phantom  Study 


Image  Type 

FBP  Reconstructed  Image 

MAP  Reconstructed  Image 

Original  Reconstructed  Image 

1.55 

1.62 

MSD  w/  Replacement  Model 

3.75 

4.05 

Table  3:  Average  Lesion  to  Background  Contrasts  for  Physical  Phantom  Study 


6.2  Dynamic  FDG-PET  Digital  Phantom  Study 

The  dynamic  digital  phantom  is  generated  by  passing  an  original  noise  free  dynamic  digital  phantom 
through  a  forward  projection  and  a  backwards  reconstruction.  The  noise  free  dynamic  digital 
phantom,  which  simulates  the  anatomic  structure  of  body  containing  heart,  normal  tissue  and 
tumors,  contains  28  frames  of  128x128  pixel  images.  To  mimic  the  clinical  physiological  process, 
all  the  TACs  assigned  to  the  digital  phantom  are  obtained  by  the  clinical  FDG-PET  studies.  The 
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FBP  Reconstructed  Image  MAP  Reconstructed  Image 


20  40  60  80  100  120  20  40  60  80  100  120 


(a)  (b) 


FBP,  MSD  w/  Replacement  Model  MAP,  MSD  w/  Replacement  Model 


(C)  (d) 


Figure  10:  Physical  phantom  study,  (a)  FBP  reconstructed  image  at  last  frame  (b)  MAP  recon¬ 
structed  image  at  last  frame  (c)  Detection  for  FBP  reconstructed  image  (d)  Detection  for  MAP 
reconstructed  image 
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measured  blood  function  of  one  patient  was  used  as  the  blood  input  function  in  the  phantom.  The 
TACs  in  the  proven  tumor  and  normal  ROI  of  the  same  patient  are  extracted  from  the  clinical 
images.  Then  based  on  the  tracer  dynamic  model  provided  by  (6),  the  parameters  are  estimated 
by  spectral  analysis  method  [13].  By  this  way,  we  can  get  low  noise  phantom  TACs  while  mimic 
true  clinic  TACs  very  well  (the  TACs  sampled  from  the  clinical  image,  even  after  averaging,  are 
with  high  noise).  The  computer  simulated  phantom  at  last  frame  was  shown  in  Figure.  6.2(a). 

When  simulating  the  forward  projection  for  the  phantom  image,  the  total  count  of  sinagram  in 
each  frame  was  scaled  to  approach  the  total  count  for  the  corresponding  frame  in  clinical  dynamic 
study.  Then  the  sinagram  data  was  added  with  Possion  noise  and  blurred  for  system  resolution. 
FBP  method  is  used  for  the  reconstruction.  Figure.  6.2(b)  shows  phantom  dynamic  image  at  the 
last  frame. 

The  estimated  lesion  and  normal  subspace  H  and  S  were  computed  using  TAC  data  sampled 
from  a  5  x  5  ROI  in  upper  lesion  in  noise  free  phantom.  Parametric  least  square  and  non-parametric 
wavelet  based  methods  were  applied  to  identify  subspaces.  Replacement  model  and  superimposed 
model  based  matched  subspace  detection  algorithm  were  employed  for  detection.  The  output  images 
are  shown  in  Figure.  6.2  (c)(d)(e)(f).  To  calculate  the  average  lesion  to  background  contrast,  mean 
lesion  and  background  values  were  achieved  by  averaging  the  5x5  lesion  and  background  ROIs 
respectively.  The  results  are  listed  in  Table  6.2. 


Image  Type 

Upper  Lesion 

Lower  Lesion 

Reconstructed  Image 

1.45 

1.05 

Replacement  Model,  LSE  subspace 

4.49 

4.59 

Replacement  Model,  wavelet  subspace 

5.01 

4.62 

Superimposed  Model,  LSE  subspace 

29.9 

17.4 

Superimposed  Model,  wavelet  subspace 

8.57 

4.50 

Table  4:  Average  Lesion  to  Background  Contrasts  for  Digital  Phantom  Study 


6.3  Clinic  Dynamic  PET-FDG  Study 

6.3.1  Protocol  and  Reconstruction  of  Dynamic  PET-FDG  Study 

The  clinical  PET-FDG  dynamic  data  was  acquired  with  a  Siemens/CTI  EC  AT  Model  953A  whole- 
body  PET  scanner.  This  device  provides  31  contiguous  transaxial  image  planes  with  an  axial  field 
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Original  Noise  Free  Phantom 


FBP  Reconstructed  Image 


Replacement  Model,  Wavelet  Subspace 
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Superimposed  Model,  Wavelet  Subspace 
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Figure  11:  Digital  phantom  study,  (a)  original  noise  free  dynamic  phantom  at  last  frame  (b)  FBP 
reconstruction  dynamic  phantom  at  the  last  frame,  (c)  detection  based  on  replacement  model,  sub¬ 
spaces  were  estimated  by  LSE.  (d)  detection  base  on  superimposed  model,  subspaces  were  estimated 
by  LSE  (e)  detection  base  on  replacement  model,  subspaces  were  estimated  using  wavelet  method, 
(f)  detection  base  on  superimposed  model,  subspaces  were  estimated  using  wavelet  method. 
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of  view  of  10.8  cm.  The  nominal  intrinsic  resolution  of  the  system  is  4  mm  in  all  3  dimensions. 
Consecutive  detector  rings  are  separated  by  tungsten  septa  in  order  to  reduce  scatter  noise.  Dy¬ 
namic  data  will  be  acquired  from  0  to  55  minute  post  injection.  The  dynamic  structure  protocol 
for  the  clinical  data  acquisition  was  listed  in  Table  1.  Two  reconstruction  methods,  namely  FBP 


Scan  Type 

Scan  Times 

Frame  Duration 

Dynamic  scan 

0-2  min 

15  sec/ frame 

Dynamic  scan 

3-5  min 

30  sec/frame 

Dynamic  scan 

6-25  min 

1  min/frame 

Dynamic  scan 

26  -  45  min 

5  min/frame 

Table  5:  Dynamic  Data  Acquisition  Protocol 


and  MAP  reconstruction,  were  used  in  our  clinical  studies. 

6.3.2  Clinical  Lung  Cancer  Study 

The  image  contains  one  confirmed  primary  tumor  and  one  confirmed  metastatic  lymph  node,  which 
is  shown  in  Figure  6.3.2  (a)  and  (b).  H  and  S  are  estimated  using  5x5  ROI  date  in  the  primary 
tumor  region.  Substituting  the  H  and  S  estimated  by  LSE  and  wavelet  methods  into  (31),  finally 
result  in  the  output  images  shown  in  Figure  6.3.2  (c)(d)(e)(f).  Notice  that  only  replacement 
model  based  methods  are  employed  in  clinical  data  study.  Because  for  clinical  data  we  have  tried, 
superimpose  model  based  methods  always  get  worse  result  comparing  to  replacement  model.  This 
problem  will  be  explained  in  the  discussion  section.  The  average  metastasis  to  background  contrasts 
are  listed  in  Table  6.3.2.  To  calculate  the  mean  contrast,  mean  metastasis  values  were  achieved 
by  averaging  the  3x3  ROI  centered  at  the  confirmed  metastasis  region,  while  mean  background 
values  were  achieved  by  averaging  whole  body’s  region(not  the  whole  image). 

6.3.3  Clinical  Breast  Cancer  Study 

The  detection  procedures  are  same  as  the  lung  cancer  case,  two  breast  cancer  images  at  different 
planes  are  tested:  one  for  the  primary  tumor,  another  for  metastatic  lymph  node.  The  original  clinic 
images  and  all  output  images  are  shown  in  Figure  6.3.3  and  Figure  6.3.3.  The  average  metastasis 
to  background  contrasts  are  listed  in  Table  6.3.2. 
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Figure  12:  Clinic  lung  cancer  study,  (a)  FBP  reconstructed  image  at  last  frame  (b)  MAP  recon¬ 
structed  image  at  last  frame  (c)  Detection  for  FBP  reconstructed  image,  subspaces  were  estimated 
by  LSE  (d)  Detection  for  MAP  reconstructed  image,  subspaces  were  estimated  by  LSE  (e)  Detec¬ 
tion  for  FBP  reconstructed  image,  subspaces  were  estimated  by  Wavelet  method  (f)  Detection  for 
MAP  reconstructed  image,  subspaces  were  estimated  by  Wavelet  method 
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Figure  13:  Clinic  breast  cancer  study,  (a)  FBP  reconstructed  image  at  last  frame  (b)  MAP  recon¬ 
structed  image  at  last  frame  (c)  Detection  for  FBP  reconstructed  image,  subspaces  were  estimated 
by  LSE  (d)  Detection  for  MAP  reconstructed  image,  subspaces  were  estimated  by  LSE  (e)  Detec¬ 
tion  for  FBP  reconstructed  image,  subspaces  were  estimated  by  Wavelet  method  (f)  Detection  for 
MAP  reconstructed  image,  subspaces  were  estimated  by  Wavelet  method 
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Figure  14:  Clinic  breast  cancer  study,  (a)  FBP  reconstructed  image  at  last  frame  (b)  MAP  recon¬ 
structed  image  at  last  frame  (c)  Detection  for  FBP  reconstructed  image,  subspaces  were  estimated 
by  LSE  (d)  Detection  for  MAP  reconstructed  image,  subspaces  were  estimated  by  LSE  (e)  Detec¬ 
tion  for  FBP  reconstructed  image,  subspaces  were  estimated  by  Wavelet  method  (f)  Detection  for 
MAP  reconstructed  image,  subspaces  were  estimated  by  Wavelet  method 


Image  Type 

Lung  Metastasis 

Breast  Metastasis 

FBP  Reconstructed  Image 

2.67 

2.03 

MAP  Reconstructed  Image 

3.99 

2.34 

FBP,  LSE  subspace 

6.84 

3.13 

FBP,  wavelet  subspace 

4.34 

2.91 

MAP,  LSE  subspace 

4.95 

5.01 

MAP,  wavelet  subspace 

6.38 

5.40 

Table  6:  Average  Lesion  to  Background  Contrasts  Using  Different  Method  for  Clinical  Lung  and 
Breast  Metastasis  Cases 


6.4  Receiver  Operating  Characteristic  (ROC)  Study 

It  is  known  that  lesion  detectability  is  well  described  by  the  ROC  curve.  The  ROC  curve  is 
determined  by  the  rates  of  true  positive  and  false  positive.  Fifty  lesion  phantom  studies  (with  250 
known  lesions)  and  fifty  normal  tissue  phantom  studies  were  conducted.  The  rates  of  true  and  false 
positive  decisions  made  from  the  three  ROI  kinetic  analyses  are  counted,  respectively,  to  form  ROC 
curves.  The  ROC  curves,  corresponding  to  the  three  kinetic  feature  analysis  methods,  are  obtained 
by  Monte  Carlo  simulations  and  shown  in  Fig.  1  (e).  The  resulting  ROC  curves  show  that  the 
ROI  based  kinetic  analyses  with  spatial  decorrelation  outperform  the  conventional  TAC  average. 
Putting  an  additional  constraint  to  the  structure  of  the  interpixel  correlation  further  lowers  the 
false  positive  rate,  because  the  correlation  structure  is  estimated  from  the  last  frame  and  applied  to 
the  rest.  Due  to  relatively  high  counts,  the  structure  estimated  in  the  last  frame  is  more  accurate 
than  that  in  the  early  frames.  The  ROC  curves,  corresponding  to  the  three  kinetic  feature  analysis 
methods,  are  obtained  by  Monte  Carlo  simulations  and  shown  in  Fig.  6.4. 

The  resulting  ROC  curves  showed  that  the  area  under  curve  (AUC)  in  the  ROI  based  kinetic 
analysis  with  spatial  decorrelation  outperform  the  conventional  TAC  average.  Putting  an  additional 
constraint  to  the  structure  of  the  interpixel  correlation  further  lowers  the  false  positive  rate,  because 
the  correlation  structure  is  estimated  from  the  last  frame  and  applied  to  the  rest.  The  reason  is 
that,  due  to  relatively  high  counts,  the  structure  estimated  in  the  last  frame  is  more  accurate  than 
that  in  the  early  frames.  Also  the  AUC  for  MAP  and  OSEM  is  larger  than  that  of  FBP,  which 
shown  that  the  iterative  reconstruction  method  outperform  the  FBP  method  not  only  in  the  image 
visualization,  but  also  in  the  detection  performance. 
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ROCstujyl 


HOC  Sim* 


Figure  15:  ROC  curves  simulation  (a)  FBP:  GLRT  (white  noise,  non- white  noise)  vs  multi-pixel 
FLRT,  where  Az (multi-pixel  GLRT)>  Az  (GLRT  with  non-white  noise)  «  Az  (GLRT  with  white 
noise),  (b)  GLRT  with  white  noise  assumption:  AZ(MAP)>  A2(0SEM)>  A2(FBP). 


6.5  Discussion 

The  phantom  and  the  clinical  data  study  show  that: 

matched  subspace  detection  result  v.s.  original  image:  MSD  results  have  higher  tumor  to  back¬ 
ground  contrast.  In  phantom  study,  the  MSD  results  can  find  lower  lesion  invisible  in  the  original 
image.  In  clinical  study,  the  MSD  results  improved  lesion  to  background  contrast  for  all  cases,  but 
for  some  case  the  high  noise  make  the  false  alarm  rate  high. 

replacement  model  v.s.  superimposed  model:  In  phantom  study,  it  is  guaranteed  that  the  pure  tumor 
TACs  can  be  expressed  explicitly  by  H  only  and  the  pure  normal  TACs  can  be  expressed  explicitly 
by  S  only.  So  the  superimposed  model,  just  as  what  we  expected,  have  higher  tumor  to  background 
contrast  in  phantom  study.  At  the  same  time  the  noise  in  the  MSD  result  with  superimposed  model 
is  also  higher  than  that  with  replacement  model.  However,  the  situation  changes  for  clinical  study. 
In  reconstructed  clinical  PET  image,  the  tumor,  not  matter  how  big  it  is,  could  be  more  or  less 
mixed  with  normal  tissue.  There  may  not  exist  such  a  pure  tumor  which  can  be  used  to  extract 
the  true  tumor  only  subspace  H.  Even  if  the  large  tumor  ROI  were  defined,  the  estimated  tumor 
subspace  is  somehow  mixed  subspace  H(JS.  So,  for  clinical  study,  the  replacement  model  provide 
better  results. 

LSE  bases  v.s.  wavelet  bases:  For  some  cases  the  LSE  bases  have  better  results,  while  for  other 
cases  wavelet  bases  provide  better  result.  However,  the  wavelet  method  always  have  one  advantage 
over  LSE  method:  LSE  needs  blood  input  function  which  is  not  always  available  or  blood  modeling 
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which  will  make  the  LS  curve  fitting  process  difficult.  On  the  contrary,  the  wavelet  method  don’t 
need  blood  input  function  all  the  time  and  is  robust  for  all  data  set  we  test. 


7  Conclusion 

We  have  obtained  encouraging  results  in  preliminary  research  of  exploiting  kinetic  features  of  the 
time  activity  curves  (TAC),  identification  subspaces  and  using  matched  subspace  filter  to  assist  in 
lesion  detection  with  FDG-PET  dynamic  images. 

Results  showed  that  the  TAC  subspaces  identified  from  the  FDG-FDG  dynamic  images  have  the 
important  characteristics  which  are  applicable  for  computer-aided  lesion  detection:  (a)  express  TAC 
curves  accurately,  (b)  distinguishable  as  lesion  and  normal  tissue  subspaces,  (c)  readily  incorporable 
to  a  matched  subspace  detector  for  lesion  detection.  By  an  application  of  the  identified  lesion  and 
normal  tissue  subspaces  to  the  matched  subspace  detector,  the  lesion-to-normal  tissue  contrast  can 
be  significantly  enhanced  relative  to  the  original  filtered  backprojection  clinical  images.  A  Monte 
Carlo  of  phantom  dynamic  study  is  conducted  for  the  detection  performance  based  a  receiver 
operating  characteristic  (ROC)  analysis. 
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A  The  condition  of  non-decreasing  slope  for  the  physiological  pro¬ 
cess  function 

Suppose  that  the  blood  input  function  Cp(t)  is  a  continuous,  differentiable  function  and  satisfies 
the  following  conditions: 

Cp(t)  =  0,  t  <  0;  Cp(t)  >  0,  t  >  0; 

dC  (t)  (37) 

c'p(t)  =  ~dt  >  °'t  <  T°'  <%(*)  =  O’1  =  X b;  c'p(t)  <  0,t  >  To, 

More  over,  we  assume  all  functions  appear  in  this  proof  are  differentiable  and  continuous. 

Let  the  physiological  process  function  be: 

/(/?,<)  =  Cp(t )  (g>  e~f3t  —  f  Cp(T)e~^t^T^dT  =  f  Cp(t  —  r)e_^rdr,  /?  >  0,  t  >  0  (38) 

Jo  Jo 

We  want  to  show  that: 

(i)  V  T  >  T0,  3/3*  s.t.  =  0 

(ii)  V  0  <  t  <  T,  >  0 

(hi)  V  /?  >  /?*,  ®^|i=T  <  0 

(iv)  V  0  <  (3  <  (3*,  a£|^|0<t<r>0 

Let  g(P,t)  =  ,  using  (38),  g(P,t)  can  be  written  as: 

g(P,  t )  =  =  Cp(t)  -  t)  (39) 

or 

g(P,  t)  =  =  Cp( 0)e-*  +  jT  C'p(r)e-^dT  =  Cp(0)e"^  +  jT  e^^(r)dr  (40) 

Proof  of  (i). 

When  fix  the  t  =  T,  g((3,T)  is  the  continious  function  of  /?.  Using  (39),  when  /3  =  0,  c/(0, T)  = 
CP(T)  >  0.  If  we  can  show  that  g(/J,  f)|t=r  <  0  for  some  /?  >  0,  then  there  must  exist  a  /?*  such 
that  s(/3*,T)  =  0. 

Let  T+  be  a  time  point  s.t.  To  <  Tq  <  T,  there  exist  a  number  0  <  M  <  oo  s.t.  Cp{t )  < 
—M,  Vt  €  [T+,T],  Also,  there  exist  a  number  0  <  N  <  oo  s.t.  Cp{t)  <  N,Vt  €  [0,Tq].  Using  the 
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later  part  of  equation  (40): 


903,  T)  =  e~0T[CP(  0)  +  [T e^C'p(r)dT\ 

Jo 

/'Tq  /*] '  ^  (J1 

=  e-^r[Cp(0)+  [  efTC'(T)dT  +  /  °  e^Cl(r)dr+  [  e0TC'{r)di 
Jo  Jt0  Jt+ 

’[Cp(0)+  /  °  epTC'{r)dT  +  fT  e^C'p{T)dT\ 

Jo  Jt+ 


<  e~PT 
<e-PT 


[Cp(0)  +  [T°e0rNdT-  C  e^M dr] 
Jo  Jt+ 


(41) 


ef3To  _  1 


e/3 T  _eHT+ 


=  e-^[Cp( 0)  +  -  M-  ~  "  6 


P 


P 


< 


o-PT 

~T 

Me-P (T-T+)  ^(0)/?  +  AT  +  M 


[Cp(0)/9e/3To+  +  NePT°  +  Me0To  -  Me131'} 


P 


M 


.  e/3(T-T0+)] 


Notice  that  T  —  7  (0  >  0,  the  term  Qs(.°)P+N+M  increases  linearly  with  (3,  while  term  e^T  To) 
increases  exponentially  with  p.  So,  for  some  P>  0,  g(P,  T)  <  0.  Now,  since  5(0,  T)  >  0  and  we  just 
proved  that  there  exist  a  value  p  such  that  g(/3 ,  T)  <  0,  there  must  be  a  point  f3*  where  g(P*,T )  =  0. 


Proof  of  (ii) 

When  fix  j3* ,  g(P*,t)  is  the  function  of  t.  By  using  (40),  it’s  easy  to  show  that  g(P*,t)  >  0  for 
t  G  (0,To].  Next,  we’ll  consider  the  case  for  t  G  (To,T).  From  (39),  taking  time-derivative  of  g(P,t) 
results  in: 


Pg(P,t)  wr.x  nPf(P,t) 

dt  p  dt 


=  C'p{t)  -  Pg(p,t) 


(42) 


We  have  assume  that  g{P*,t)\t=T  =  0.  If  g(P*,t)  is  not  always  positive  in  (T0,T),  suppose  at  time 
point  i!  G  (7b,  T),  g(P* ,  tr)  <  0.  Then  there  must  exist  another  time  point  t"  G  (t',T)  such  that: 
d9<'gt’t‘>  |t=t"  >  0,  C'p(t")  <  0  and  g(P*,  t")  =  0,  which  contradict  (42).  So  fort  G  (0,T0),  g(P*,t)  >  0. 
Combining  two  above  cases,  V0  <  t  <  T,  —  1^*’^  >  0 


Proof  of  (iii) 

In  order  to  prove  (iii),  we  first  need  to  prove  that  \p=g*  <0,  t  G  (0,T].  Continue  from 
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(40) 


=  -tCp{0)e~f*  -  j\t  -  T)C>p(T)e-W-^dr  (43) 

This  equation  shows  that  <  0  is  hold  for  any  f3>  0  and  t  G  (0,  To].  Proceed  from  (43): 

=  — crp(0)e-*  +  ptCp( 0)e"*  -  jT  C'p{T)e-^dr  +  pj\t~  r)C'(r)e-^dr 

=  -Cp(O)e~0t  -  f  C'(t )e~P('t~T'> dr  +  /?[tCp(0)e-#  +  [\t  -  r)Cl(r)e-«*-T)dr] 

•to  ./o  (44) 

_WM)  _  a2/(/M) 

at  p  dm 


To  prove  9  gpgft  \ p=p*  <  0  for  te  (0,  T),  the  prove  by  contradiction  is  used.  Suppose  there  exist 
a  solution  to  \p=p*,t=t0  =  0,  to  G  (To,T).  Without  loss  the  generality,  we  can  also  assume 

the  to,  if  it  exists,  is  the  minimum  one.  So,  9  g^g'J  <  0  for  /?  =  /?*,  t  G  (To,  to).  Using  equation(44) 


d3/(/M)  , 

ata/^t 


g/(/g,*) 

dt 


\p-p-,t=t0  <  0 


(45) 


Notice  that  to  6  (0,T),  by  using  (ii),  the  ”<”  sign  is  hold.  Remind  that  <  0  for  t  €  (0,To]. 

As  t  increase  from  To,  due  to  our  assumption,  to  is  the  first  time  point  that  the  function 
becomes  zeros.  This  implies  that  ^Jpgif  >  Q\p=p*,t=t0,  contradict  with  (45).  So,  there  does  not 
exist  to  such  that  9  \p=p*,t=t0  =  0  for  t  G  (0,T).  Therefor 

<0’  *G(0>T)  (46) 

Next,  we’ll  show  9  g^§f^  \p=p*,  t=T  <  0.  Using  equation  (43), 

d2Q{J~i%=0>,t=T  =  -TCP( 0)e-P'T  -  J\t  -  r)C'p(T)e-^T-^dr 

=  ~T[CP(  0)e-^T+  fT  C'JT)e~fi'{‘T-TUT]+  [T  rCL(T)e~ ^T~r^  dr  (47) 

Jo  Jo 

=  -Tg(P*,T)+  [T  rC'p{r)e-^T^ dr 
Jo 
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Notice  that  g(/3*,T)  —  0,  so: 


d2f(/3,t) 

d(3dt 


\p=p\t=T  =  [  TC'p(T)e  0*{T  TW 
Jo 

—  [  TC'(T)e~0{'T~T'>dT+  f  TC’p(T)e~l3*^T~T^dT 

Jo  Jt0 

<  [T°  ToC'p(T)e-0^T-^dT  +  (T  nC'p(T)e-0t(?-TUT 
Jo  J  To 

=  T0  F  C'p(r)e-0*V-^ dr 
Jo 

=  T0[g(f3*,T)  -  Cp( O)e~0'T}  =  -T0Cp(0)e~^T  <  0 


So,  d qm/'i-\ p\  t=T  <  o.  Combine  this  with  (46),  we  get: 


d2f(M 


1/3=/?*  <  0)  *€(0,T] 


Next,  we’ll  prove  V/3  >  /?*,  \t=T  <  0  by  using  (49)  and  prove  by  contradiction.  If  (iii) 

is  not  hold,  then  when  (3  increases  from  (3*,  assume  the  function  \t=T  cross  the  zero  first 

time  when  (3  =  (3*+ .  Notice  that,  same  as  /?*,  /3*+  also  satisfies  (ii)  and  (49).  This  implies, 
I /3=p*+,  t=T  <  0-  But  this  contradicts  with  the  assumption  that  /3*+  is  the  first  value  which 
makes  \t-T  cross  the  zero.  There  for: 


V  0  >  /?*, 


9f(P,  t ) 


\—T  <  0 


Proof  of  (iv) 

Use  the  same  technique  as  proof  of  (iii),  we  can  prove  V  0  </?</?*,  2^|t=r  >  0.  So,  we  only 
need  to  prove  that 

VO  </?</?*,  ^^|o<t<r>0  (51) 

The  proof  by  contradict  will  be  employed  again.  Assume  (51)  is  not  hold,  because  is 

continuous  function,  there  must  exist  a  (3*~,  0  <  (3*~  <  (3*  and  a  T~,  0  <  T~  <  T  such  that 
d^dt’t'>  1/3=/?*-,  t=T~  =  0-  Then  (3*~  and  T~  are  a  pair  of  f3  and  t  values  which  satisfy  (i),  (ii)  and 
(iii).  To  make  to  clear,  rewrite  (iii)  for  /?*“  and  T~  case: 

v /?>/?*',  <  0  (52) 
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Notice  that  /3*  >  f3*  ,  T  <  T,  above  equation  implies  that  p=p*,  t=T~ 

obviously  contradict  with  (ii).  So,  equation  (51)  is  hold.  Finally,  result  in: 


<  0  which  is 


VO  <P</3*, 


dt 


|o <t<T  >  o 


(53) 
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Abstract 


Residual  FDG  activity  in  normal  tissues,  such  as  blood  vessels  and  the  liver,  as  well  as  the  spilled 
in  background  activity  can  impair  the  detection  of  small  or  modestly  tracer-avid  tumors  in  FDG- 
PET  cancer  imaging.  In  order  to  maximize  tumor  visualization  but  minimize  background  and 
artifacts,  an  efficient  new  method  is  adapted  from  the  constrained  temporal  filtering  processing 
widely  used  for  signal  detection  and  extended  to  the  application  of  dynamic  FDG-PET  processing. 
Comparing  with  the  well-known  Patlak  analysis  and  spectral  analysis  (SA),  the  main  advantage 
of  the  proposed  method  is  that  it  can  objectively  remove  the  partial  volume  effect  superimposed 
onto  the  TAC  as  well  as  the  residual  blood  activities  and  result  in  the  pixel-by-pixel  estimations  of 
the  infux  constant  as  the  output  of  the  filter.  Since  the  constrained  temporal  filter  is  designed  to 
preserve  the  power  of  tumor  signal,  therefore,  it  is  likely  to  offer  a  desirable  noise  canceling  result, 
while  make  no  or  minimum  distortion  on  tumor  signatures.  In  contrast,  Patlak  analysis  may  suffer 
in  their  ability  to  classify  tumor  and  normal  structures  when  the  tumor  is  severely  interfered  by 
background  activity  and  present  approximately  similar  time  activity  curve  as  that  of  normal  tissues 
at  the  end  of  the  post-FDG  injection  period.  For  comparison,  analytical  expressions  for  SNR  and 
CNR  of  filtered  images  are  derived,  and  the  results  of  our  digital  phantom  study  suggest  that 
the  proposed  method  has  the  best  performance,  comparing  with  the  Patlak  and  SA.  The  mice 
implanted  with  cancer  cells  are  imaged  with  dynamic  FDG-PET  on  the  6th  day  and  the  12th  day 
after  implantation,  respectively.  The  proposed  method  is  applied  to  process  the  earlier  acquired 
images  and  the  resulting  findings  are  confirmed  with  the  later  acquisitions.  The  results  show  that 
the  new  method  can  enhance  the  tumor-to-background  ratio  at  an  early  stage  and  is  promising  for 
improving  lesion  detectability. 


1  Introduction 


Positron  emission  tomography  (PET)  with  [18F]2-fluoro-2-deoxy-D-glucose(FDG)has  been  widely 
used  for  non-invasive  diagnosis  of  cancer.  The  diagnostic  advantage  of  FDG-PET  is  important 
because  alterations  in  the  metabolism  of  cells  are  often  evidenced  in  disease  states  before  struc¬ 
tural  changes  can  be  determined  using  other  medical  imaging  technologies,  such  as  computerized 
tomography  (CT)  and  magnetic  resonance  imaging  (MRI).  While  PET  imaging  provides  a  higher 
sensitivity  (95)  for  detecting  lymph  node  metastasis  than  alternative  modalities  [1],  recent  study 
shows  nearly  a  20  false-negative  rate[2].  In  addition  to  the  relatively  limited  spatial  resolution  of 
PET,  the  residual  blood-born  FDG  activity  may  be  contributing  to  this  finding.  Note  that  most 
dynamic  FDG-PET  tumor  detection  methods  are  based  on  one  fact  that  malignancies  can  be  dis¬ 
tinguished  from  normal  tissues  on  the  basis  of  biologically  determined  radiotracer  accumulation 
or  loss  rate.  Because  malignant  tumors  are  metabolically  active  and  FDG-avid  on  PET  imaging, 
they  metabolize  glucose  at  a  much  higher  rate  than  do  most  normal  tissues  [3].  However,  even  after 
60  minutes  post-FDG  injection,  there  is  substantial  blood-pool  FDG  activity  remaining,  possibly 
impairing  the  ability  to  detect  small  tumor  near  blood  vessel  or  near  tissues  with  high  normal  FDG 
activity.  Besides  the  effect  of  residual  FDG  activity,  the  partial  volume  effect  (PVE)  also  imparis 
the  spatial  resolution  of  PET.  PVE  is  contributed  from  the  surrounding  background  activities,  the 
boundary  pixels  of  a  malignant  area  represent  a  mixture  of  tumor  and  normal  tissues,  instead  of 
the  pure  tumor  tissues.  Therefore,  PVE  dramatically  lowers  the  tumor-to-background  contrast  of 
small  structures  that  are  smaller  than  two  times  the  full-width  at  half-maximum  (FWHM)of  the 
tomograph  [4]. 

Developing  methodologies  for  maximally  removing  the  residual  blood  FDG  activity  and  cor¬ 
recting  the  PVE  has  always  been  a  critical  issue  in  various  FDG-PET  quantitative  studies.  For 
instance,  the  non-parametric  techniques  simply  decompose  the  observed  FDG-PET  data  to  extract 
the  desired  signal  subspace  using  singular  value  decomposition  or  eigen-decomposition  method. 
After  kinetic  modeling  using  FDG  is  well  established,  parametric  methods  attract  more  atten¬ 
tion.  The  widely  used  Patlak  analysis  estimates  the  macroparameter  by  determining  the  slope  of  a 
transformed  tracer  uptake  curve,  which  is  based  on  three-compartment  three-rate  constants  model 
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(3K)  [5] .  This  method  is  extremely  computationally  efficient  and  has  good  performance  on  residual 
FDG  activity  cancellation  [2] .  However,  the  assumption  underlying  the  Patlak  method  is  that  FDG 
is  irreversibly  trapped  in  the  system(/c4  =  0)  and  this  assumption  only  be  valid  for  most  malignant 
tumors  with  progresive  FDG  accumulation  in  tissue  over  time.  Therefore,  Patlak  analysis  may  be 
vulnerable  when  the  tumor  is  severely  interfered  by  background  activity  because  of  PVE.  Spectral 
analysis  is  another  conventional  parametric  approach,  which  enables  the  determination  of  pharma¬ 
cokinetic  parameters  with  relatively  few  model  assumptions  [6].  This  method  is  limited  by  a  lack  of 
knowledge  of  parameters  in  the  model  and  the  fitting  procedure  is  also  a  concern,  which  is  based  on 
non-negative  lease  squares  (NNLS)  algorithm,  not  a  linear  operator.  Other  previous  attempts  to 
produce  high  quality  PET  images  considered  the  use  of  various  filters  for  FDG-PET  images  either 
in  the  spatial  or  frequency  domains [7], [8], [9]. 

This  paper  proposes  a  new  approach  by  using  a  constrained  temporal  filter.  The  emphasis 
of  this  method  is  to  preserve  the  FDG  features  attributed  to  tumor  tissues  and  eliminate  the 
undesired  contributions  of  the  residual  FDG  activity  at  the  same  time.  In  fact,  by  using  two 
imposed  constraints,  we  can  objectively  remove  the  partial  volume  effect  superimposed  onto  the 
time-activity  curves  (TACs)  as  well  as  the  residual  blood  activities,  while  making  no  distortion  on 
the  desired  malignancy  features  in  time  domain.  Moreover,  because  this  constrained  temporal  filter 
is  designed  to  minimize  the  energy  of  the  output  signal,  our  approach  can  also  suppress  the  PVE 
and  statistical  noise,  besides  removal  of  the  structured  interferences.  Simulation  results  in  section 
IV  show  that  the  proposed  approach  outperforms  Patlak  and  spectral  analysis. 

The  paper  is  structured  as  follows.  Sections  II  introduce  the  FDG  model  and  briefly  review 
Patlak  and  spectral  analysis.  Sections  III  presents  the  ...  Section  IV  validates  ...  through  phantom 
and  animal  studies. 

Notation:  All  boldface  letters  indicate  vectors  (lower  case)  or  matrices  (upper  case). 

2  The  Fluoro-Deoxy-Glucose  (FDG)  Kinetic  Model 

2.1  A  4-k  Compartmental  Model  for  the  FDG  Tracer  Kinetics 

The  FDG  kinetic  model  consists  of  three  compartments  which  describes  the  tracer  distribution 
in  the  homogeneous  tissue,  as  shown  in  Figure  1.  The  first  compartment  is  considered  as  the 
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blood  pool  region  and  represents  the  concentration  of  FDG  in  the  plasma.  The  second  and  third 
compartments  contain  the  concentrations  of  FDG  in  tissue  and  phosphorylated  FDG  in  tissue 
(FDG-6-P),  respectively.  The  rate  constants  in  Figure  1  are  defined  as  follows: 


K\  =  the  transport  of  FDG  in  plasma  to  tissue, 
fc2  =  the  transport  of  FDG  in  tissue  to  plasma, 
&3  =  the  phosphorylation  of  FDG  to  FDG-6-P, 
k4  =  the  dephosphorylation  of  FDG-6-P  to  FDG 


The  rate  of  change  in  the  concentrations  of  FDG  and  FDG-6-P  in  a  homogeneous  ROI  can  be 
described  by  the  following  differential  equations 


dCe(t) 

dt 

dCmjt) 

dt 


K\Cp(t)  —  (k,2  +  ks)Ce(t)  +  k4Cm(t) , 


face(t)  -  k4Cm(t), 


(1) 


where  Ce(t)  represents  the  concentration  of  FDG  in  the  tissue  at  time  t,  Cm(t )  represents  the 
concentration  of  FDG-6-P  (and  all  metabolites  derived  from  FDG-6-P)  in  the  same  region,  and 
Cp(t)  represents  the  FDG  concentration  in  the  arterial  plasma.  Solutions  of  (1)  [10]  given  that  the 
initial  conditions  Ce(0)  =  Cm(0)  =  0  show  that  the  total  radioactivity  for  a  homogeneous  tissue, 
C(t),  is  the  sum  of  the  free  [18F]FDG  concentration  plus  the  concentration  of  metabolites,  i.e., 


where 


C(t)  =  Ce(t)  +  Cm(t) 

=  [(*»  +  h  -  Pi)e +  {(32  -  h  -  k4)e-^}  ®  Cp{t) 

P2  ~  PI 

=  [. M ie~^  +  M2e-^‘]  g>  Cp(t), 
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(k2  +  ^3  +  k4)  —  V  (fc2  +  k^  +  k4)2 
{k2  +  &3  +  k4)  +  yj  (fc2  +  k%  +  k4)^ 
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Plasma  Cell  Tissue  Cell 

Figure  1:  4-k  three  compartment  model 


and  <g)  denotes  the  convolution  operator. 

The  radioactivity  measured  by  the  PET  scanner  (under  noise-free  assumption),  including  the 
radioactivity  in  the  cerebral  blood,  can  be  represented  as 

<W*)  =  Vbcp(t)  +  C(t)  =  VbCp(t)  +  [Mie"*‘  +  M2e-0*t]  ®  Cp(t ),  (4) 

where  CtQtei\(t)  is  the  total  radioactivity  measured  by  the  scanner  at  time  t;  Vb  denotes  the  vascular 
space  in  tissue  and  FDG  concentration  in  whole  blood  is  assumed  to  be  VbCp(t)  at  time  t  [11]. 

2.2  A  3-k  Compartmental  Model  for  the  FDG  Tracer  Kinetics 

The  3-K  compartment  FDG  model  is  proposed  by  Sokoloff  el  al.  in  1977  [10].  This  model  assumes 
that  [18F]FDG-6-P  is  irreversibly  trapped  in  tissue  for  the  duration  of  experiment  and  it  is  described 
by 

SS  =  KiCp(t)  -  Ik,  +  k,)CJt), 

=  hC,(t),  (5) 

The  total  radioactivity,  including  the  radioactivity  in  the  cerebral  blood,  is  represented  as 

CtoUt)  =  Vbcp(t)  +  C(t)  =  VbCp(t )  +  [Ml  +  M2e-^]  ®  Cp(t),  (6) 

3  CONVENTIONAL  METHODS  OF  DYNAMIC  PET  TUMOR 
DETECTION 

Before  presenting  the  proposed  constrained  temporal  filtering  approach,  it  is  worthwhile  to  present 
a  brief  overview  of  Paltlak  and  spectral  analysis  in  general.  We  will  compare  the  performance  of 
these  two  methods  with  that  of  the  proposed  approach  in  Section  V. 
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3.1  Patlak  analysis 


Patlak  analysis  assumes  3-K  compartment  model  for  the  FDG  tracer  kinetics  (that  is,  k\  is  assumed 
to  be  zero). 

The  solution  to  the  differential  equations  (3)  is 

rt 


CtotaiW  =  Ki  f  Cp{r)dT  +  vCp(t) 

Jo 


(7) 


where  v  is  related  to  the  effective  distribution  volume  of  FDG. 

By  dividing  through  by  Cp(t),  the  Patlak  method  employed  to  determine  the  macroparameter 
is  given  by 

Ctotai (t)  _  K  fo  Cp(r)dT 

cp(t )  -Ki  Cp(t)  +  v  w 


Based  on  above  equations,  Patlak  approach  is  implemented  by  an  pixel-by-pixel  analysis  of 
^  Cp(t) ' ^  versus  )^' ’  w^ich  gives  the  estimation  of  Ki.  That  ATj-slope  image  is  expected  to  be 

a  one  parameter  approach  to  improve  image  contrast  over  static  imaging,  since  the  tumor  would 
have  a  higher  influx  and  retention  of  FDG  over  most  other  normal  tissues  as  a  function  of  time.  In 
particular,  earlier  work  in  [2]  shows  that  Patlak  analysis  has  good  performance  on  residual  FDG 
activity  cancellation. 

However,  Patlak  analysis  may  be  vulnerable  because  the  assumption  underlying  the  Patlak 
method  is  that  FDG  is  irreversibly  trapped  in  the  system  (fc4=0).  Even  this  assumption  is  valid 
for  most  cases  because  most  malignant  tumors  have  progressive  FDG  accumulatio  in  tissue  over 
time,  it  is  failed  sometimes  when  the  tumor  is  severely  interfered  by  background  activity  because 
of  PVE.  Figure  5(b)  shows  the  time-activity  curves  (TACs)  extracted  in  the  selected  ROIs  from  a 
mouse  study  at  an  early  stage  of  tumor  growth.  As  we  can  see  the  TAC  of  the  tumor  tissue  has  an 
approximately  flat  slope  at  the  later  frames.  If  we  apply  Patlak  analysis  to  compute  the  A^-slope 
image,  it  fails  to  classify  the  tumor  and  normal  tissues. 


3.2  Spectral  analysis 


In  the  spectral  analysis  method,  the  tissue  concentration  of  FDG  activity,  CUlUil(t),  is  modeled  as  a 
linear  combination  of  basis  functions,  each  defined  as  a  exponential  function  convolved  with  Cp(t ): 
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t  >  0 


J 

<W«)  =  aie~Pit  ®  cp{t), 

j=i 
j 

t~° 
j= i 

where  J  is  the  maximum  number  of  basis  functions  allowed  in  the  model  and  values  of  f3  are  cho¬ 
sen  to  cover  the  spectrum  of  expected  kinetic  behavior.  A  spectrum  of  (3  values  are  predefined  and 
spaced  logarithmically  on  a  given  interval  [A,l],  where  A  is  the  decay  constant  of  the  radioisotope. 
The  unknown  coefficients,  ay ,  can  be  estimated  by  the  linear  least  square  recursive  curve  fitting 
technique  with  non-negative  constrains.  This  technique  can  be  applied  at  a  pixel  level  enabling 
calculation  of  parametric  images  independent  of  an  assumed  compartmental  structure  [11].  There¬ 
fore,  it  is  a  general  modeling  approach.  However,  due  to  a  lack  of  knowledge  on  the  range  of  /?, 
covering  the  spectrum  of  all  possible  kinetic  behavior  often  results  in  estimation  of  a  large  number 
of  unknown  coefficients  ays.  The  accuracy  of  the  estimation  is  then  degraded.  Furthermore,  the 
linearity  of  the  spectral  analysis  method  may  break  down  if  (a)  one  or  more  of  the  measurement 
vectors  contains  negative  values,  or  (b)  the  true  number  of  non-zero  coefficients  cannot  be  resolved, 
either  because  of  measurement  noise  or  insufficient  degrees  of  freedom.  [12] 

4  CONSTRAINED  TEMPORAL  FILTER 

4.1  Optimization  Criteria 

Based  on  3-K  compartment  model  and  taking  into  account  the  partial  volume  effect,  we  rewrite 
the  equation  (7)  as 

Ctotalit )  =  Ki  f  Cp(r)dr  +  VbCp(t)  +  Cv(t)  +  N{t)  (10) 

J  o 

where  Cv(t)  is  the  sum  of  ®  C*  (t)]  and  partial  volume  effect  (PVE)  contributed  from  the 

surrounding  pixels  (voxels)  ;  N(t)  is  the  random  noise  which  is  uncorrelated  with  Cp(t)  and  Cv(t). 
Let  vector  cp  =  [Cp(ti),  Cp{t 2), . . . ,  Cp(tm)}T  be  the  sampled  plasma  concentration,  where  tj  is  the 
time  instant  that  frame  j  is  acquired,  and  similarly  define  the  vectors,  c totaU  cv  and  n.  Equation 
(10)  then  be  expressed  in  a  vector  form, 
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C total  —  -E^Sp  4“  -f-  C.v  "F  n  —  A)®p  d"  6 


(11) 


where 

ftl  pt2  ftm 

s  P  =  [  Cp(r)dT,  /  Cp(r)dT, . . . ,  /  Cp(t  )dr  (12) 

Jo  Jo  Jo 

e  =  Vbcp  +  cv  +  n  (13) 

In  order  to  cancel  the  interference  e  defined  in  equation  (13),  pervious  work  has  developed 
many  ways  to  find  the  weight  vector  w  of  an  optimum  filter  subject  to  well-defined  optimality 
criteria.  Although  the  details  differ  in  different  applications,  the  main  assumptions  and  processing 
algorithms  are  essentially  the  same.  In  particular,  an  important  assumption  is  that  the  interfering 
signals  are  uncorrelated  with  the  desired  signal.  However,  in  our  practice,  the  signal  (KiSp)  and 
the  interference  (e)  may  be  either  partially  or  completely  correlated.  Moreover,  [13]  shows  how 
the  interference  rejection  becomes  more  severe  with  increasing  correlation  between  the  desired 
signal  and  the  interferences.  Therefore,  in  this  work,  we  introduce  a  new  optimum  criterion  called 
optimum  interference  plus  noise  rejecter  [14].  By  using  this  criterion,  all  the  interferences,  whether 
correlated  or  uncorrelated  with  the  desired  signal,  are  considered  undesired  and  the  weight  vector  w 
is  chosen  to  pass  the  /Q  undistorted  at  the  filter  output  while  maximally  rejecting  the  contribution 
of  the  residual  FDG  activity,  PVE  and  random  noise.  Multiplying  both  sides  of  equation  (11)  with 
w,  the  output  of  the  optimum  filter  is  given  by 


WTC  total  =  KiVfTSp  +  VbWTCp  +  VfTCv  +  WTn 


(14) 


=  KiWTsp  +  wTe 

where  wTe  is  the  undesired  contribution  of  the  residual  FDG  activity,  PVE  and  noise. 

Then  we  propose  two  constrains  to  ensure  unity  gain  in  the  desired  signal  direction  and  complete 
removal  of  the  interference, 


wrsp  =  1  (15) 

wTe  =  0  (16) 
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However,  because  interference  e  includes  unknown  random  noise  n  and  background  FDG  con¬ 
centration  bs,  it  is  usually  difficult  to  determine  the  exact  format  of  the  e.  In  practice,  we  define 
another  interference  variable  esuf,  as  esuf,  =  cp  or  esu/,  =  [cp,bs],  which  is  depended  on  whether  we 
can  successfully  extract  the  features  of  the  interference  caused  by  bg.  Then,  we  can  rewrite  (16)  as 

~wTesub  =  0  (17) 

and  the  optimization  objection  function  can  be  expressed  as 

min  £[(wre)(wTe)T]  =  min  wTRew(18) 

W^Sp  =  1  wrs  p  =  1 

wT©j mb  —  0  wTesub  =  0 

where  Re=E(eeT)  is  the  covariance  matrix  of  the  interference  and  it  is  free  of  the  desired  signals. 
Substituting  the  constraints  (15)  and  (17)  into  (14),  the  output  of  the  optimum  filter  is 

Ki  =  •wTctotai  =  Ki  +  wTct)  +  wTn  «  Kt  (19) 

The  approximation  in  (19)  holds  because  the  interference  e,  including  random  noise,  should  be 
minimized  by  the  filter  w  resulted  from  Equation  (18). 

4.2  Constrained  Filter 

The  addition  of  the  constrains  to  adaptive  processing  was  first  described  by  Frost  [15].  Frost 
presented  a  method  of  adapting  the  weights  in  a  filter  while  holding  the  filter  response  fixed  in  a 
certain  waveform  or  time  function  or  at  certain  frequencies  [16].  It  means  that  if  we  have  a  mxl 
constraint  vector  h,  defined  in  prior  by  a  filter  designer,  no  matter  how  the  filter  weights  Wi  are 
changed,  we  can  always  keep 

wTh  i  =  fi  (20) 

where  fi  is  a  given  constant  and  w  is  a  mxl  weight  vector. 

Equation  (12)  illustrates  that  the  constrained  filter  will  preserve  all  the  energy  of  hi  if  fi  equals 

X 

to  unity  and  eliminate  the  energy  of  hi  if  fi  equals  to  zero.  Suppose  there  are  K  such  constraint 
equations  altogether.  Then  we  may  combine  these  constraints  in  one  equation  by  defining  a  mx  K 
constraint  matrix, 

H  =  [h!,h2, . . . ,  h.K]  (21) 
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and  by  defining  f  to  be  a  A-dimensional  vector  of  the  coefficients  fi 


f=[/i,/2,...,/ic]T  (22) 

The  constrains  may  then  by  written  as 

wtH  =  fT  (23) 

If  we  let  the  interference  e(13)  be  the  input  signal  of  the  constrained  filter,  then  the  output  signal 
is, 

e  =  WTe  (24) 

The  average  power  of  the  constrained  filter  output  signal  is 

P  =  wTf?(eer)w  =  wTRew  (25) 

Hence  the  optimization  problem  is  to  minimize  wTReW  subject  to  the  constraint  (23).  Comparing 
with  (18),  that  is  exactly  the  optimization  criterion  what  we  started  in  the  previous  part.  By  using 
LMS  criteria  [18],  we  have 

w  =  Re-1H(HrIU-1H)-1f  (26) 

which  is  the  weight  vector  of  our  constrained  temporal  filter,  satisfying  the  optimization  criterion 
(18).  The  derivation  of  (26)  is  given  in  Appendix  A. 


In  the  following  discussion,  we  will  use  the  above  derivation  to  find  the  the  constrained  temporal 
filter  in  order  to  enhance  the  tumor-to-background  contrast  in  dynamic  PET  imaging.  Suppose  the 
image  size  is  Lx L  pixels  at  a  given  frame.  Combining  all  images  of  m  frames  together,  we  form  a 
new  mxL2  matrix, 


X  —  [C-total  j  C total  ^total  :  C total  ] 


(27) 


where  c totail  is  the  TAC  in  vector  form  at  the  ith  pixel  as  illustrated  in  (11). 

Even  though  our  goal  is  to  minimize  wTReW  (18)  (25),  it  is  actually  difficult  to  determine  the 
value  of  interference  e  and  its  covariance  matrix  Re  since  they  include  the  unknown  random  noise 
n  and  background  FDG  concentration  bs.  In  practice,  we  take  and  alternative  way  to  achieve  that 
goal. 


9 


Because  the  tumors  under  our  study  are  assumed  being  small,  in  the  TAC  observations  c total1 , 
the  residual  FDG  activity,  PVE  and  noise  should  be  dominant,  comparing  with  the  FDG  activity 
of  tumor  tissue.  Therefore,  we  can  use  samples  of  c totail  at  all  different  pixels  to  approximate  the 
covariance  matrix  of  esub, 


® sub  : 


’  X  —  [ctotal  i  C total  j  -  •  •  j  C total  )  d total  ] 


(28) 


hence 

I C,  «  Re  *  XXT  (29) 

and  using  the  two  constrains  stated  in  (15)  and  (16),  we  construct  the  constraint  matrix, 


H  —  [Sp,  ©sub] 


(30) 


and 


/  [1)  0])  if  &sub  —  Cp  ! 

\  (1))  if  esub  =  [Cp,  bg]  . 


f  —  [1 5  0] ,  if  Gsub  —  Cp 

or  f=  [1,0,0],  if  esub  =  [cp,bg] 


(32) 


Substituting  (30),  (33)and  (30)  into  (26), we  have  the  optimum  constrained  temporal  filter  for 
dynamic  PET,  which  satisfied  (15),  (16)  and  (18) 

4.3  Extension  to  4-K  model 

We  can  also  adapt  the  proposed  optimum  criterion  to  the  4-K  compartment  model.  Using  the  same 
definitions  as  in  (6),  the  tissue  concentration  of  FDG  activity  is  rewritten  as: 

c totality  =  Mi  [  e_/?10-r)Cp(r)dr  +  Vbc p(t)  +  c„(t)  +  n (t)  (34) 

Jo 

and  the  vector  form  is  given  by: 

Ctotai  =  MiSp  +  He p  +  Cv  +  n  =  M\sp  +  e  (35) 
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where 


ft\  pt2  ptm 

sp  =  [  /  e-/3i(ti-r)c p(r)d,T,  /  e-/3i(i2-'r)Cp(r)dr, . . . ,  /  e~^tm~T^c p(r)dr  (36) 

Jo  Jo  Jo 

Then,  we  can  follow  the  same  derivations  as  those  of  3-K  compartment  model. 

In  general,  performance  of  tumor  detection  is  degraded  with  respect  to  that  of  the  3-K  model. 
This  is  expected,  since  the  4-K  model  has  additional  parameter  with  respect  the  3-K  model.  How¬ 
ever,  how  to  determine  that  additional  parameter  Pi  of  4K  model  is  a  difficult  problem.  Previous 
work  has  proposed  many  methods  to  resolve  this  problem,  like  classical  nonlinear  least  squares 
(NLS),  weighted  integration  method  (WIM)  and  generalized  linear  least  squares  (GLLS)  method. 
For  most  cases,  those  techniques  work  well  when  the  tumor  location  is  known.  However,  if  the 
tumor  is  unknown,  we  have  to  go  through  every  pixel  in  the  image  to  estimate  the  parameters. 
It  is  time  consuming  and  the  computational  requirement  is  high.  Moreover,  how  to  set  the  initial 
value  is  also  a  strong  concern.  Therefore,  in  the  proposed  approach,  we  take  a  different  way  to 
implement  the  parameter  estimation. 

It  is  shown  in  [20]that  cp(f)  satisfies  two  constrains  in  the  noise  free  case:  first,  c p(t)  >  0 
for  t  >  0;  second,  there  exists  a  constant  Tk  such  that  >  0  when  t  <  Tk  and  <  0 

when  t  >  Tk .  With  these  constrains,  it  has  been  proved  that  for  a  given  time  T,  T  >  Tk,  there 
exists  a  threshold  /?*  of  the  macroparameter  in  dynamic  FDG-PET  for  each  individual  patient, 
such  that  the  physiological  process  function  specified  by  p{  is  horizontal  (e.g.  9cpW®e  =  q) 
at  a  certain  time  t  =  T  >  Tk-  Further  more,  when  Tk  <  t,  the  physiological  functions  with 
P  <  P\  will  increase  monotonically,  or  otherwise  decrease  when  Pi  >  p\.  We  also  notice  that 
the  the  physiological  process  function  for  tumor  cell  will  go  up  along  with  time;  physiological 
process  function  of  inflammation  cell  will  be  roughly  fiat  after  a  long  time;  while  the  normal  cell’s 
physiological  process  function  will  go  down.  Therefore,  wehn  decomposing  a  TAC  into  a  linear 
combination  of  many  physiological  process  functions,  we  hypothesize  that  the  physiological  functions 
with  Pi  smaller  than  p\  are  associated  to  malignancy  tissue,  but  the  physiological  functions  with  Pi 
larger  than  P\  are  related  to  normal  tissue. 

Based  on  above  finding,  we  can  first  determine  the  value  of  /?*,  which  is  the  first  step  of  our 
parameter  estimation  method.  Then,  we  choose  several  predefined  P\  values  by  spacing  the  interval 
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[0,  /?*]  uniformly  or  logarithmically.  In  the  next  step,  we  form  a  constrained  filter  bank,  in  which 
each  filter  is  corresponding  to  one  predefined  value  of  0 1  and  generated  by  the  proposed  approach 
of  this  work.  Suppose  there  are  N  branches  in  the  filter  bank,  after  using  this  constrained  filter 
bank  to  process  the  dynamic  PET  images,  the  output  vector  of  a  certain  pixel  in  the  PET  image 
is  given  by 

r=  [r0,n,...,r-jv-i]r  (37) 

n  =  w,Tc  total  =  M1wtTsp  +  w  iTe  (38) 

=  Si  +  Zi  (39) 

where  s;  is  the  desired  signal  of  the  ith  branch  in  the  constrained  filter  bank  and  Z{  is  the 
interference  and  noise  of  branch  i. 

Because  z  =  (zo, —  l]T  is  not  white,  the  covariance  matrix  Rz  =  zzT  is  not  necessarily 
a  scaled  identity  matrix.  Therefore,  we  pass  the  output  vector  r  through  a  noise  whitening  filter, 
the  resulting  signal  is  given  by: 

r/  =  s/  +  z/  =  sRz-1/2  +  zRz-1/2  (40) 

where  r/  =  [ro',  rR,  . . . ,  rjy_ \l)T . 

Finally,  we  use  Maximal  Ratio  Combining  (MRC)  technique  to  combine  the  output  of  all 
branches  simultaneously.  At  this  stage,  each  of  the  branch  output  is  weighted  with  a  gain  fac¬ 
tor  proportional  to  its  own  SNR  and  all  of  the  output  are  combined  to  give  the  resulting  output 

image,  which  has  the  best  SNR  value.  The  output  of  the  MRC  scheme  can  be  written  as  follows: 

N- 1 

y=Yl 9iri'  (41) 

i= 0 

where  gi  is  the  combining  factor  of  branch  i  and  the  derivation  of  gt  is  given  in  Appendix  B. 

The  constrained  filter  bank  and  maximal  ratio  combining  scheme  is  shown  in  Figure.  2. 

How  to  extract  the  features  of  the  interference  caused  by  background  FDG  concentration  bg 
is  another  strong  concern  of  the  proposed  technology.  The  simplest  way  to  resolve  this  problem  is 
that  we  do  not  use  any  background  information  when  we  form  the  esub,  ie.  ,  esub  =  cp.  However,  we 
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Figure  2:  Constrained  temporal  filter  bank 


still  want  to  eliminate  the  background  interference  as  much  as  possible,  since  the  background  FDG 
concentration  will  definitely  degrade  the  image  quality  of  the  detection  result.  In  our  approach, 
we  draw  ROIs  in  the  dynamic  PET  image  and  obtain  the  bs  from  the  ROI  of  the  background 
area.  Then  we  form  the  esub  as  esub  —  [cp,  bs](33)  and  generate  the  optimum  weight  vector  w 
for  the  whole  image.  The  best  way  to  form  the  isub  is  to  acquire  the  features  of  hg  adaptive.  For 
instance,  we  can  use  computed  tomography  (CT)  to  segment  the  whole  PET  image  and  distinguish 
several  different  background  area.  Then  we  draw  ROIs  in  these  background  area  and  extract  several 
difficult  features  bj,  b^,. .  ,,b”.  The  background  interference  is  expressed  by 

esub  =  [cP>kg]  (42) 

then  we  use  elsub  to  generate  the  optimum  weight  vector  w1  for  the  ith  segmentation  of  the  PET 
image. 

4.4  Decision  Criterion 

The  problem  of  detecting  small  areas  of  images  which  differ  in  some  statistical  sense  from  their 
immediate  surroundings  is  of  considerable  interest  in  a  number  of  image  processing  applications. 
This  section  addresses  a  power  detection  algorithm  which  is  the  locally  optimal  detector  for  random 
objects  in  Gaussian  backgrounds.  For  the  3-K  compartment  model,  the  two  hypotheses  of  our 
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problem  are  given  by  (14): 


Ho  :  Xi  =  ni—  wTe 

(43) 

Hi  : 

Xi  =  Si  +  nt  =  K,wTsp  +  wre 

(44) 

where  Xi  is  the  a  certain  pixel  in  the  ith  segmentation  of  the  PET  images,  s*  is  the  equivalent 
desired  signal  and  rij  is  the  equivalent  noise.  We  assume  that  rq  is  Gaussian  process  with  local 
mean  mi  and  power  spectral  density  of. 

Because  signal-plus-noise  density  is  unknown  due  to  the  uncertainty  of  tumor  location,  likelihood 
ratio  test  can  not  be  applied  in  our  case.  However,  since  we  can  measure  the  first  and  second 
moment  statistics  of  equivalent  noise  rii  from  the  images,  power  detection  scheme  can  provide  great 
probability  of  detection  for  a  given  false  alarm  rate.  This  power  detector  consists  of  a  half-wave 
linear  device  with  the  transfer  characteristic 


f  Xi,  when  Xi>  0  ;  , 

^  =  \0,  when  x;  <  0  .  (45) 

Then  by  taking  the  power  of  y{  and  applying  a  power  detector,  the  decision  criterion  becomes 


"if  Zi  >  then  decide  object  present ” 


(46) 


with  the  threshold  /(a*)  is  a  function  of  false  alarm  rate  a* 


Ui  =  Q(- 


CTi 


where  Zi  is  the  power  of  yi  and  the  Q-function  is  defined  by 


(47) 


1  /*°° 

Q(X)  =  —J I  e-‘/2df  (48) 

The  derivation  of  Equation  (48)  is  given  in  Appendix  C. 

Following  the  same  procedure,  we  can  derive  the  decision  criterion  for  the  4-K  compartment 
model. 


5  EXPERIMENT  RESULTS  AND  DISCUSSION 

5.1  Digital  phantom  Study 

In  our  digital  phantom  study,  we  use  the  simulated  dynamic  PET  images,  containing  a  uniformly 
distributed  ellipsoid  background  and  9  small  tumors.  Among  them,  the  three  tumors  in  the  upper 
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row  are  crosses  of  3  x  2  pixels,  containing  5  pixels  in  total.  The  middle  and  bottom  three  tumors  are 
squares  of  2  x  2  pixels  and  3x3  pixels,  respectively.  The  area  of  each  pixel  is  2.7 mm  and  tumor- 
to-background  ratio  is  equal  to  1.5.  The  time-activity  curves  of  the  tumor  and  normal  tissue, 
which  are  used  to  generate  this  original  phantom  image,  are  extracted  from  the  clinic  patient 
data.Three  different  methods  have  been  applied  to  process  the  phantom  image:  Patlak  analysis, 
spectral  analysis  (SA)  and  the  constrained  temporal  filter.  Visual  results  Figure. 3(a)  shows  the 
original  phantom  image  (true  image),  which  contains  nine  artificial  tumors  in  a  uniform  backgroud. 
Two  noisy  and  noise-free  Filtered  Backprojection  (FBP)  reconstructed  phantom  images  are  showed 
in  Figure. 3(b)  and  Figure. 3(c).  Both  of  figures  show  the  last  frame  of  the  reconstructed  phantom 
images,  which  are  believed  to  have  the  highest  tumor-to-background  contrast  among  al  dynamic 
images.  By  comparing  the  Figure. 3(a)  and  Figure. 3(b),  we  can  clearly  see  the  distortion  caused  by 
partial  volume  effect.  Figure. 3(d),  Figure. 3(e)  and  Figure. 3(f)  are  the  results  of  the  noisy  image  in 
Figure. 3(c)  processed  by  Patlak,  SA  and  the  constrained  temporal  filter,  respectively.  Figure. 3(g) 
shows  the  TACs  drawn  from  the  tumor  and  background  in  the  original  phantom  image.  Figure. 3 
(h)is  the  TACs  drawn  from  the  boundary  of  the  tumor  in  the  noise-free  FBP  reconstructed  phantom 
image,  which  actually  is  the  mixture  of  the  two  TACs  showed  in  the  Figure. 3(g)  because  of  PVE. 

5.1.1  Quantitative  results 

For  performance  evaluation,  we  measured  the  SNR  and  CNR  in  the  output  images  resulted  from 
the  above  three  methods.  The  definitions  of  SNR  and  CNR  are  described  next. 

(1)  Singal-to-Noise  Ratio 

We  use  the  same  definitions  of  SNR  and  CNR  given  in  [19].  Let  Vjkd  be  the  value  of  the  (j,  k)th 
pixel  in  the  desired  ROI(DROI).  Vjkd  consists  of  a  deterministic  value  plus  statistical  noise;  the 
desired  deterministic  value  if  then  the  mean  E(Vjkd),  and  the  strength  of  the  noise  is  the  standard 
deviation  [var(Vjkd)^2  ■  Signal-to-noise  ratio  of  the  desired  feature  ( SNRd )  is  defined  as 

SNRd  = 

where  E(Vjkd)  and  [var{Vjkd)^2  in  a  homogeneous  region  can  be  estimated  using  the  sample  mean 
and  variance. 


15 


Figure  3:  Results  of  digital  phantom  study,  (a)  Original  phantom  image,  (b)  FBP  reconstructed 
phantom  image  (noise  free),  (c)  FBP  reconstructed  phantom  image  (with  noise),  (d)  Noisy  image 
in  Figure. 3(c)  processed  by  Patlak  analysis,  (e)  Image  in  Figure. 3(c)  processed  by  SA,  (f)  Image  in 
Figure. 3(c)  processed  by  constrained  temporal  filter. 
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TABLE  I 

PERFORMANCE  OF  PATLAK  ,  SA  AND  CONST.  FILTER 


SA 

Patlak 

Const.  Filter  1 

SNR 

CNR 

SNR 

CNR 

SNR 

CNR 

Tumor  (1,1) 

1.157 

0.833 

2.238 

1.472 

1.434 

1.064 

Tumor  (1,2) 

1.672 

2.009 

1.113 

0.954 

1.385 

1.403 

Tumor  (1,3) 

0.945 

1.213 

2.184 

2.090 

3.298 

3.326 

Average  of 
above  3  tumors 

1.258 

1.352 

1.845 

1.505 

2.039 

1.931 

Tumor  (2,1) 

1.322 

1.468 

1.270 

0.755 

2.136 

1.789 

Tumor  (2,2) 

1.169 

1.414 

2.084 

1.848 

3.172 

3.026 

Tumor  (2,3) 

0.500 

0.254 

1.700 

1.037 

3.609 

3.487 

Average  of 
above  3  tumors 

0.997 

1.045 

1.685 

1.214 

2.972 

2.767 

Tumor  (3,1) 

0.618 

0.240 

1.359 

1.068 

1.633 

1.472 

Tumor  (3,2) 

1.185 

1.418 

2.054 

1.788 

3.188 

3.001 

Tumor  (3,3) 

1.715 

1.960 

3.086 

2.700 

2.854 

2.724 

Average  of 
above  3  tumors 

1.173 

1.206 

2.166 

1.852 

2.559 

2.399 

(2)  Contrast-to-Noise  Ratio 

In  addition  to  SNR^,  the  CNR  between  the  desired  signal  (tumor  tissue)  and  undesired  signal 
(interference)  is  usually  important  for  image  interpretation.  The  CNR  between  the  desired  signal 
and  undesired  signal  is  defined  as 


CNRdu  — 


E(Vjkd)  -  E(Vjku) 


(50) 


yJo.5[var{Vjkd)  +  var(Vjku)] 
where  Vjku  is  the  value  of  the  (j,  k)th  pixel  in  the  undesired  ROI  (UROI),  and  E(Vjku )  and  var(Vjku ) 
are  the  mean  and  the  variance  of  pixel  values  in  UROI,  respectively. 

By  using  (49)  and  (50),  Table  I  compares  the  quantitative  results  of  Patlak  analysis,  spectral 
analysis  (SA)  and  the  constrained  temporal  filter.  As  Fig. 3(a)  shows,  the  nine  tumors  in  original 
phantom  data  are  arranged  as  a  3  x  3  matrix.  The  index  (i,j)  of  that  matrix  are  used  in  Table  I 
to  indicate  the  location  of  tumors. 

Both  image  results  and  numerical  results  suggest  that  the  constrained  temporal  filter  has  the 
best  performance.  Comparing  with  the  other  two  methods,  the  main  advantage  of  our  approach 
is  that  it  can  objectively  remove  the  PVE  superimposed  onto  the  TAC  and  the  residual  blood 
activities  in  the  pixel-by-pixel  estimations  of  M\ . 
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5.2  Physical  phantom  Study 


Figure  4:  Results  of  physical  phantom  study. 

(a)  MAP  reconstructed  phantom  image  (b)  Output  image  of  processing  Figure.4(a)by  the  constraint 
temporal  filter,  (c)  Power  detector  result  of  Figure.4(b), 


(a)  (b)  (c) 


(a)  MAP  reconstructed  phantom  image  (b)  Output  image  of  processing  Figure. 4(a) by  the  constraint 
temporal  filter,  (c)  Power  detector  result  of  Figure. 4(b), 

5.3  Animal  Study 

Dynamic  PET  imaging  of  FDG  (200  fid)  was  performed  on  a  MicroPET  R4  system  (Concorde 
MicroSystems,  Inc).  Female  athymic  nude  mice  with  human  breast  cancer  cell  MDA-MB-435 
inoculated  into  both  sides  of  mammary  fat  pads  under  the  arms.  Six  days  after  inoculation,  thirty- 
five  dynamic  data  frames  (6x1  sec,  4x3  sec,  10  x  30  sec,  5  x  60  sec,  and  10  x  300  sec)  were  acquired 
for  1  hour  after  intravenous  injection.  Images  were  reconstructed  with  the  OSEM  algorithm,  as 
shown  in  Figure. 5(a).  The  corresponding  time-activity  curves  of  the  two  tumors,  normal  tissues 
and  heart  have  already  been  shown  in  Figure. 5(c).  Figure. 5(e)  is  the  resulted  image  of  Figure. 5(a) 
after  passing  through  the  constraint  temporal  filter. 

In  Figure. 5(b)  we  show  the  OSEM  reconstructed  image  acquired  twelve  days  after  inocula¬ 
tion.  Figure. 5(d)  demonstrates  the  corresponding  TACs  and  Figure. 5(f)  is  the  resulted  image  of 
Figure. 5(b)  after  passing  through  the  constraint  temporal  filter.  Comparing  Figure. 5(a)  and  Fig¬ 
ure. 5(b),  we  can  see  that  the  tumors  grow  bigger  and  the  slope  of  tumor’s  TAC  becomes  larger 
after  six  more  days.  Based  on  detection  result  of  Figure. 5(h),  we  confirm  that  the  findings  in  the 
filtered  output  image  shown  in  Figure. 5(g)  are  correct. 
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From  the  results  of  the  mouse  study,  it  appears  that  using  constrained  temporal  filter  is  an 
optimal  approach  to  eliminate  the  residual  blood  FDG  activity  especially  for  the  stressing  cases 
with  difficulties  to  tell  the  difference  between  the  slopes  of  tumor  and  normal  tissues.  This  approach 
also  offers  a  desirable  result  in  eliminating  partial  volume  effect  (PVE)  and  image  noise  artifacts. 

5.4  Clinic  Study 

5.4.1  Protocol  of  Dynamic  PET-FDG  Study 

The  clinical  PET-FDG  dynamic  data  was  acquired  with  a  Siemens/CTI  ECAT  Model  953A  whole- 
body  PET  scanner.  This  device  provides  31  contiguous  transaxial  image  planes  with  an  axial  field 
of  view  of  10.8  cm.  The  nominal  intrinsic  resolution  of  the  system  is  4  mm  in  all  3  dimensions. 
Consecutive  detector  rings  are  separated  by  tungsten  septa  in  order  to  reduce  scatter  noise.  Dy¬ 
namic  data  will  be  acquired  from  0  to  55  minute  post  injection. The  dynamic  structure  protocol  for 
the  clinical  data  acquisition  was  listed  in  Table  2. 


Scan  Type 

Scan  Times 

Frame  Duration 

Dynamic  scan 

0  -  2  min 

15  sec/frame 

Dynamic  scan 

3-5  min 

30  sec/frame 

Dynamic  scan 

6-25  min 

1  min/frame 

Dynamic  scan 

26  -  45  min 

5  min/frame 

Table  1:  TABLE  II:  Dynamic  Data  Acquisition  Protocol 


5.4.2  Clinical  Breast  Cancer  Study 

Two  breast  cancer  images  at  different  planes  are  tested:  one  for  the  primary  tumor,  another  for 
metastatic  lymph  node. 

5.4.3  Clinical  Lung  Cancer  Study 

The  detection  procedures  are  same  as  the  breast  cancer  case.  The  image  contains  one  confirmed 
primary  tumor  and  one  confirmed  metastatic  lymph  node 


19 


0. Normal 

Tissue 

20- 

30  - 

\ffr 

40  - 

/  |  \ 

KLeft 

1  Right 

Tumor 

M- 

Heart  Tumor 

Normal 

Tissue 

**  1 

•  -  X 

so  Left  > 

Right 

Tumor  Heart  Tumor 

4  ,  ;  ; 

_ . _ . _ o 

10  20  30.  .  40  50  50  10  »  30  -  40  SO  GO 

(a)  (b) 


Figure  5:  Results  of  animal  study,  (a)  OSEM  Reconstructed  FDG-PET  image  (six  days  after 
inoculation)  at  the  last  frame  with  the  highest  tuflfi>r-to-backgournd  ratio,  (b)  OSEM  Reconstructed 
FDG-PET  image  (twelve  days  after  inoculation)  at  the  last  frame  with  the  highest  tumor-to- 
backgournd  ratio,  (c)  TACs  drawn  from  Figure.5(a),  (d)  TACs  drawn  from  Figure. 5(b),  (e)  Output 
image  of  processing  Figure.4(a)by  the  constraint  temporal  filter,  (f)  Output  image  of  processing 
Figure. 4(b)by  the  constraint  temporal  filter,  (g)  Power  detector  result  of  Figure.5(e),  (h)  Power 
detector  result  of  Figure.5(f). 


6  CONCLUSION 


Enhancement  of  tumor-to-background  contrast  in  dynamic  PET  images  may  be  obscured  by  resid¬ 
ual  FDG-associated  activity  in  normal  tissues  and  blood  as  well  as  by  PVE.  Even  though  Patlak 
method  and  spectral  analysis  can  improve  tumor  visualization  when  tumors  are  large,  it  fails  some¬ 
times  when  tumor  and  normal  tissues  have  the  almost  same  slope  of  TAC  due  to  the  PVE.  A  new 
approach  of  using  the  constrained  temporal  filter  seems  to  be  superior  to  Patlak  and  SA  methods 
as  shown  in  our  digital  phantom  study  and  the  animal  study  of  early  detection  with  mouse. 
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7  APPENDIX 

7.1  A 

The  problem  of  optimization  subject  to  a  constraint  may  be  handled  with  the  Lagrange  multiplier 
method.  We  define  the  performance  index 


I  (w)  =  wJ  Rew  +  A(wi  H  —  r ) 


(51) 


where  A  is  A'-dimcnsional  column  vector  of  Lagrange  multiplier.  We  may  determine  the  optimum 
weight  vector  w  by  setting  the  derivatives  of  /( w)  with  respect  to  Wi  equal  to  zero, 


Then,  we  obtain 


where 


-  =  0, 1  <  i  <  m  and  w  =  \wi,W2, . . . ,  wm}T 

OWi 


Vw/(w)  =  Rew  +  HA  =  0 


d  d 


d 


L<9u>i’  dw2’  ’  dw„ 


The  optimum  weight  vector  is  then 


w  =  -Re  XHA 


(52) 

(53) 

(54) 

(55) 


To  evaluate  w,  we  have  to  know  A.  However,  since  w  must  satisfy  the  constraint,  the  following 
relation  holds: 

Htw  =  -HTRe_1HA  =  f  (56) 

Hence, 

A  =  -(HTRe-1H)_1f  (57) 

and  Equation  (55)  gives 

w  =  Re-1H(HrRe-1H)-1f  (58) 

7.2  B 


The  combined  output  of  the  MRC  scheme  can  be  written  as  follows: 

N- 1  N- 1  N- 1 

y=Yl 9iri' =  XI 9iSi' +  giZi' 


(59) 


i=0  i=0 


i= 0 
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where  gi  is  the  combining  factor  of  branch  i  and  zj  of  the  different  branches  are  uncorrelated.  The 
SNR  of  the  combined  output  signal  is  given  by 

r  =  (60) 

E*=o  9i2°i 

where  cq2  is  the  power  spectral  density  of  ztt  and  Ef  is  the  energy  of  ry.  The  combining  gain 
factors  should  be  chosen  in  order  to  maximize  the  above  SNR.  By  using  Schwarz  inequality  we  can 
write  the  above  resulting  SNR  as 

r  <  (Eiln1  gjWKES ElMl  y  El 

-  a.2a.2  ~  2^  a  2 
2^i=0  9i  ai  i=0  * 

Equality  in  Schwarz  inequality  is  obtained  if  and  only  if  combining  gain  factors  are  chosen  as 

9i  =  k - 2  (62) 

where  k  is  an  arbitrary  positive  constant.  Thus,  the  maximum  value  of  the  output  SNR  is  given  by 

N- 1  rp  i  N-l 

r=E^=Er  '<  (63> 

i= o  1  i= 0 

Consequently,  the  SNR  of  the  resulting  signal  is  equal  to  the  sum  of  the  SNRs  in  each  branch  for 
maximal  ratio  combining. 

7.3  C 

For  the  desired  signal  absent  hypothesis  Ho,  the  probability  density  function  of  x,  is  given  by 

1  (xj-mj)2 

Pxi(xi)  =  f—  e  (64) 

V27T<7, 

The  probability  distribution  function  of  the  output  of  the  half-wave  linear  device  can  be  expressed 
in  terms  of  the  input  borbability  density  function  as 

P(vt  <Vi)=0  for  Vi<  0  (65) 

fVi 

or  P(yt  <  yi)  =  P{xi  <  0)  +  /  pXi(xi)dxi  for  yi>  0  (66) 

Jo 

By  differentiating  Equation  (66)with  respect  to  yi,  we  get 


PyiiVi )  =  Pfa  <  0 )6{yi)+pXi{yi)U{yi) 
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where  U(yi)  is  the  unit  step  function.  Since  Zi  —  (yi)2,  yi>  0 


pM  =  Q(-r)s(^)  +  rrr--6  (68) 

The  false  alarm  rate  a*  is  given  by 

Oii=  pZi(Zi)dZi  =  Q( - )  (69) 
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